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DEVELOPMENT OF A SENSITIVITY ANALYSIS TECHNIQUE 
FOR MULTILOOP FLIGHT CONTROL SYSTEMS" 


ABSTRACT 

Traditionally, the process of developing the aerodynamic model 
of an aircraft has been done independently from the process of 
designing the aircraft control system. The result has been that, 
very often, additional wind-tunnel and flight tests are performed to 
refine the model beyond what is required for control system 
design. This separation between the process of refining the 
aerodynamic model and the design of the control system has led to 
inefficient use of resources. 

This report presents the development and application of a 
sensitivity analysis technique for multiloop flight control 
systems. This analysis yields very useful information on the 
sensitivity of the relative-stability criteria of the control 
system, with variations or uncertainties in the system and 
controller elements. This technique could be used to integrate the 
modelling and controller design processes, through the determination 
of the parameters of the aircraft and controller models that have 
the largest effect on the relative stability of the closed-loop 
system. 
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The sensitivity analysis technique developed is based on the 
computation of the singular values and singular-value gradients of a 
feedback-control system. The method is applicable to single- 
input/single-output as well as multiloop continuous-control 
systems. The use of the singular values to obtain a relative 
stability criterion for sampled-data systems is also explored in 
this report, with promising results. A possible way to extend the 
sensitivity analysis method for sampled-data systems is also 
suggested in this report but has not been implemented. 

The sensitivity analysis technique was applied to a continuous 
yaw/roll damper stability augmentation system of a typical business 
jet, and the results show that the analysis is very useful in 
determining the system elements which have the largest effect on the 
relative stability of the closed-loop system. 

As a secondary product of the research reported here, the 
relative stability criteria based on the concept of singular values 
were explored. This application to a single- loop yaw damper system 
and a multiloop yaw/roll damper system of a business jet was 
demonstrated. The results showed that the stability margins 
obtained through the singular- value analysis are conservative when 
compared to the margins obtained from a Bode analysis. A way to 
reduce the conservativeness of the analysis is also discussed. 
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1 . INTRODUCTION 


1.1 BACKGROUND 

During the development of a flight control system, a series of 
refinements of the mathematical model of the aircraft dynamics are 
used, each corresponding to a different stage of the vehicle 
development. Often, additional wind-tunnel and flight tests are 
performed to refine the model beyond what is required for control 
system design. This separation between the process of refining the 
model and the design of the control system has led to inefficient 
use of resources. 

One approach to the problem of designing a controller using a 
model which may not accurately represent the aircraft dymamics is to 
design a controller which is insensitive to model uncertainties or 
variations. This is the general approach known as "robust 
controller design." A disadvantage of using this approach is that 
the final controller design may be unnecessarily "robust," and this 
could translate into greater complexity and cost or reduced 
performance. 

An alternative approach to the problem would be to integrate 
the process of model determination with the process of controller 
design, via a sensitivity analysis. Since the modelling process 
requires resources to be spent on theoretical estimates and on wind- 
tunnel and flight tests, a knowledge of how the controller 
performance changes with changes in the elements of the model would 
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provide a way to direct those resources to the accurate 
determination of only the important controller and aircraft system 
parameters. This selected allocation of resources would be expected 
to increase the efficiency of the overall system analysis and 
synthesis problem. 

An integrated procedure. Figure 1.1, would start with the 
development of a preliminary model that would be based on 
theoretical calculations and available wind-tunnel data or flight 
test data. With the use of this model, a controller is designed 
which can be either analog or digital and can be designed using 
classical or optimal design procedures. Regardless of the details 
of the design, all of the designs will have some degree of 
robustness or insensitivity to variations in the model parameters. 

A sensitivity analysis is performed to determine the effect of model 
parameter variations on controller performance. This will allow an 
ordering of the relative importance of accurate estimates of the 
individual model parameters. Those model parameters that produce 
large changes in the closed-loop system performance for small 
changes in the parameter values may warrant more accurate 
determination if those changes in performance endanger the stability 
of the closed-loop system or degrade the handling qualities of the 
aircraft. With the aid of this sensitivity information, specific 
tests could be designed to determine accurately the important 
parameters. If, in the testing, the changes in the values of the 
parameters were found to be significant, the controller could be 
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Figure 1.1: Integrated Procedure for Controller Design 

and Model Determination 
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redesigned. This procedure offers the potential for reducing the 
amount of wind-tunnel testing, flight testing, and data analysis 
required. 

The present report covers the research done under NASA 
Cooperative Agreement NCC 2-293, at the Flight Research Laboratory 
of the University of Kansas, for the development of a multiloop- 
system sensitivity analysis method that could be used to integrate 
the modelling and controller design processes. 

As a result of the evaluation of the state of the art of 
sensitivity analysis methods for multiloop systems, a method based 
on the singular values was selected. Also developed during this 
work was a method of calculating the gradients of singular values 
with respect to model parameters. The singular values and singular- 
value gradients were combined to form a sensitivity analysis 
technique. The use of the singular values as a multi loop- system 
analysis tool is demonstrated by their application to some simple 
third-order systems, a single-input/single-output yaw damper, and a 
multiloop yaw/roll damper system for a typical business jet. The 
use of the sensitivity analysis method developed is also 
demonstrated by its application to a simple, third-order system and 
to the multi loop yaw/roll damper system for a business jet. 
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1.2 REPORT OVERVIEW 


The remainder of this report is divided into six chapters and 
four appendices. Chapter Two presents a summary of the different 
sensitivity analysis techniques, which are discussed in greater 
detail in Reference [29]. Chapter Three introduces the concept of 
singular values, and explores the characteristics of the singular- 
value analysis as compared to traditional frequency response 
methods. This is done through the analysis of several single- 
input/single-output (SISO) and multiloop (MIMO) systems. Chapter 
Four presents the mathematical development of the singular-value 
gradients and discusses the sensitivity analysis method developed. 
The application of the sensitivity analysis method is also explored 
in Chapter Four through the use of examples. Chapter Five 
summarizes the results obtained and presents some recommendations 
for future research, and Chapter Six contains the literature 
references. Appendix A contains a user's manual for the computer 
program used to calculate the singular- value gradients. Appendix B 
contains the aerodynamic and geometric data for the typical business 
jet used in the examples, and Appendix C contains a suggested 
approach to the development of singular-value gradients for sampled- 
data systems. 
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2. REVIEW OF SENSITIVITY ANALYSIS TECHNIQUES: 
ALTERNATIVE METHODS 


The key element in the integration of the modelling process 
with the controller design processs is the sensitivity analysis of 
the control system. The development of the sensitivity analysis 
involves several fundamental questions? i.e., its applicability to a 
multi loop system, the type of mathematical model used (state 
variable form or transfer function form) , and the controller 
performance measure (such as absolute stability, relative stability, 
and handling qualities). 

The different techniques for sensitivity analysis were grouped 
as six different approaches, according to the performance measure 
used, their applicability to multi loop systems, and their use in 
discrete or continous systems. As it should be expected, not all 
the approaches are totally different from each other? and some of 
them serve as building blocks for others. These techniques have 
been discussed in Reference [29] , and only a brief summary is 
presented in Table 2.1. 

The classical sensitivity approach (Approach No. 1) refers to 
the "percentage change" sensitivity function for SISO systems? i.e., 
the ratio of the percentage change of system transfer function to 
percentage change of process transfer function or parameter. The 
eigenvalue sensitivity approach (No. 2) refers to the change in the 
eigenvalues of a system due to variations or inaccuracies in the 
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Table 2.1: Characterization of 




Applicable 

System 

Sensitivity 

Approach 

Models 

Measure 

1 . 

Classical 

• Linear 

Ratio of % change 


Sensitivity 

• Continuous 

of system TF to 



• SISO 

% change of process 
TF or parameter 

2. 

Eigenvalue 

• Linear 

Gradient of 


Sensitivity 

• MIMO, SISO 

eigenvalues with 



• Continuous 

model parameters 

3. 

Performance 

• Linear 

Performance Index 


Index 

• MIMO, SISO 

Variation, 6l 


Sensitivity 

• Continuous 



4. Trajectory • Linear a) State Vector 

Sensitivity • MIMO, SISO Gradient 

• Continuous b) Comparison of 

trajectory error 
between closed- 
loop and 
equivalent 
open- loop system 

5. Gain • Linear Variation of 

Sensitivity • Digital optimal gains 

• MIMO, SISO with parameter 

variation 


Notes: TF = Transfer Function PM = Phase Margin 

SV = Singular Value GM = Gain Margin 


Different Sensitivity Analysis Approaches 


Comments 

• SISO Systems 

• Limited applicability to complex 
MIMO problems 

• In the MIMO case there is no 3, 4, 5, 6 

clear relation to PM and GM. 

• The relevance of the variation is measured 7,8,9,10/11 

by comparing the closed-loop system with an 

equivalent open-loop system* 

• Practical applicability is very limited: 

a) The designer is usually interested in the 
optimal values and sensitivities of the 
components of "I” rather than of "I" itself. 

b) This approach is not applicable for the case 
where "I" is independent of _X_( t ) • 

c) This approach is not applicable in cases 
where the end condition jUt-j ) is not free. 

• The relevance of the variation would be 12,13,14,8,9 

determined from the comparison of open-loop 

and closed-loop sensitivities. 

• No indication of what parameter variation 
is more important. 

• Requires testing of a matrix function of the 
return difference for positive definiteness. 

This was regarded (Ref. 16) as too complex 
to be useful as a design tool. 

• Effect of sampling time T included 1 5 

• Possibility of coupling this approach with 
approach 6 could be evaluated. 


References 

1,2 




Table 2.1: Characterization of the Different Sensitivity Analysis Approaches (continued) 


Approach 

Applicable 

System 

Models 

Sensitivity 

Measure 

Comments 

References 

6. a-Plots, 

a-Gradients, 
and Universal 
Phase-Gain Margin 
Evaluation 
Diagram 

• Linear 

• Digital, 
Continuous 

• MIMO, SISO 

Gradient of SV 
with parameter 
variations + 
a-Plot + PM 
and GM Diagram 

• Relative importance of each parameter 
variation is determined using the SV gradients. 

• The relevance of the parameter variation 
change is determined using the a-plot and the 
PM and GM diagram. 

• Software for a-plots exists. 

16,17,18,19, 
20,21 ,22,23, 
24,25,26,27 


• Direct relationship with PM and GM. 

• Results are conservative. 

• Current trend in control theory. 

• New developments appearing. 

• Research work being done at 
Langley Research Center. 


Notes: TF = Transfer Function 

SV = Singular Value 


PM = Phase Margin 
GM = Gain Margin 



parameters of the system. 'The performance index sensitivity 
function (Approach No. 3) examines the changes in the performance 
index of an optimal closed-loop control system when a parameter 
vector deviates from its nominal value, and compares them with the 
changes in the performance index of the equivalent optimal open-loop 
control system. In general terms, the trajectory sensitivity 
(Approach No. 4) refers to the change in the state vector components 
with variation of system parameters. The two broad categories under 
this approach are the state-vector gradients and the comparison 
sensitivity. The first one allows for the explicit computation of 
derivatives of the root mean square value of particular state-vector 
components with respect to given parameters. The second one 
compares the state- vector trajectory error between closed-loop and 
equivalent open-loop control systems. The gain sensitivity approach 
(Approach No. 5) deals with the combined effects of sample rate and 
large parameter variations on the feedback gains of discrete flight 
control systems. Approach No. 6 utilizes the concepts of singular 
values and their related multiloop gain and phase margins (sigma 
plots) and combines them with the variation of the singular values 
with changes in system parameters (sigma gradients). 

A goal of our research was to use a performance measure for the 
control system that could be easily interpreted by control system 
designers. The gain and phase margins criterion for SISO systems 
represents the type of performance measure for which the control 
system designer has already developed a good "feeling." This is due 
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to the vast experience with the design of SISO systems that has been 
gained through the years. There fore , during the evaluation of the 
different techniques, emphasis was placed on the possible connection 
of the sensitivity analysis information and the relative stability 
of the system. For the remainder of this report, we will not be 
concerned with any type of performance measure other than the 
relative stability of the closed-loop system. 

The method developed uses the singular values and singular- 
value gradients of the return difference matrix of a feedback 
control system. The use of the singular values was selected because 
of their direct applicability to multi loop (as well as SISO) 
systems, their applicability to discrete and continuous systems, the 
availability of efficient software to evaluate the singular values, 
and the direct relationship of the singular values with multiloop 
phase and gain margins. 

Mathematical expressions for the gradients of the singular 
values with respect to system and controller parameters were 
developed; and, together with the sigma plots and the concept of 
multiloop relative stability margins, they were used to obtain a 
multiloop sensitivity analysis technique. This technique is 
discussed in Chapter Four, along with the development of the 
singular-value gradients. The fundamental concepts of singular 
values are introduced in Chapter Three. 
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3. SINGULAR-VALUE ANALYSIS 


3.1 FUNDAMENTAL CONCEPTS 

In this section, some fundamental concepts about singular 
values are introduced: their mathematical definition, their 

interpretation from the point of view of the control system analyst, 
and their relationship to multiloop gain and phase margins, 

A general multiloop feedback control system is represented in 
Figure 3.1, where G(s) and H(s) are matrices of transfer functions 
called transfer matrices. The Laplace transform of the state- 
variable vector, X + , can be written as 

X = [(I + GH) -1 G]R (3.1 ) 

where R is the Laplace transform of the command-input vector, and I 
is the identity matrix. The control-input vector, U, can be written 
as 

U = [(I + HG ) ” 1 ] R (3.2) 

where, in both equations, the argument (s) has been omitted for 
convenience in writing. The system in Figure 3.1 can be represented 
by either Equation (3.1) or Equation (3.2). The matrices (I + GH) 
and (I + HG) are called the return difference matrices of the system 
for the loop broken at the plant output and input, respectively. 


t 


Vectors are indicated in boldface type. 
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Defining as "absolutely stable" a system which will produce a 
bounded output as response to a bounded input, it can be seen that 
instability is related to the condition that (I + GH) or (I + HG) 
are singular. If (I + HG) or (I + GH) are "nearly singular," the 
system will be sensitive to modelling errors [22], since in such a 
case, a model error could make the true system unstable even though 
the model of the system is stable. Therefore, it is also important 
to consider a measure of relative stability of the system. This 
measure can be given by the "near singularity" of (I + HG) or 
(I + GH). 

There are several ways to determine if a matrix A is 
singular. One such way is the determinant of the matrix; i.e., if 
det [ A] = 0 (3.3) 

then matrix A is singular. However, the determinant is not a good 
indicator of near singularity [22]. The closeness of a matrix to 
singularity is related to the "size" of that matrix. Therefore, a 
different kind of test is needed which provides a more reliable 
measure of the "size" of the matrix. In the control system context, 
a physically meaningful measure of the "size" of a system described 
by a transfer matrix ought to bear some relation to the input-output 
gain ratio of the system [17], since such a transfer matrix 
transforms input signals into output signals. 

Consider a matrix, A, which performs a linear transformation of 
a vector, r, into a vector, u (Figure 3.2). We are interested in 
defining the measure of the "size" of A that is related to the 
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effect of the transformation of r into u. To eliminate the linear 
effect of the size of r on its image, we are interested in the value 
of II Aril /Hr II , the factor by which the size of r is stretched under 
the transformation by A. The maximum and minimum stretching factors 
can be taken as definitions of the "size" of the matrix, A: 


II All 

max 


II Aril 

II r II max 


0(A) 


(3.4) 


II All . = 

min 


II Ar II 

llr II min 


a(A) 


(3.5) 


where a (A) and a (A) are the maximum and minimum singular values of A 

[17, 20]. The maximum and minimum singular values of a matrix, A, 

are given by the positive square root of the maximum and minimum 
* * 

eigenvalues of A A, respectively, where A is the conjugate 
transpose of A (the transpose of A if A is a real matrix). 

For example if 

13 5 11 

_ * 

A = then A A = 

24 11 25 

★ 

and the eigenvalues of A A are 

= 29.85 ; \ 2 = 0.15 

and, therefore, the singular values of A are 
a = positive /29.85 = 5.46 
a_ = positive /0. 15 = 0.39 

The calculation of the singular values is generally not carried out 
by first computing the eigenvalues of A*A. An efficient algorithm 
for the numerical computation of the singular value decomposition 
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( SVD ) is generally used [35], which yields the singular values and 
their corresponding orthonormal eigenvectors. 

If a(A) = 0, then A is singular. Therefore, the minimum 
singular value of A, oiA) , can be interpreted [16] as the distance 
between the matrix and the nearest singular matrix. Also note that, 
by definition, a(A) _> 0. Applying these concepts to the system of 
Figure 3.1, a condition for relative stability is 

£(I + ) > 0 (3.6) 

where G^ — HG. 

In order to examine the system gain and phase margins, a square 
diagonal matrix, L, is introduced at the plant input, as shown in 
Figure 3.3(a); a convenient representation of the same system, when 
analyzed with the loop broken at the plant input, is shown in Figure 
3.3(b). Matrix L has the special form 

j<f> n 

L = DIAG[k e ] (3.7) 

n 

At nominal conditions, k n = 1 and <j> n = 0 for all n; and therefore L 
is the identity matrix. By selecting k n / 1 or <l> n / 0, it is 
possible to introduce a gain or phase change in loop n; this can be 
done for all loops simultaneously. The stability conditions are 

£( I + G^ ) > 0 (3.8) 

and 

Oil + G^ L) > 0 (3.9) 
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for the nominal and perturbed systems, respectively. It can be 
shown that the stability of the perturbed system is guaranteed if 

odf 1 - I) < £(I + G 1 ) (3.10) 

assuming L -1 exists [20, 22]. Lehtomaki, et al., [22] show that if 
inequality (3.10). holds, and 

£(I + G i) JL a (3.11) 

for some constant, a _<_ 1 , then simultaneously in each loop of the 
feedback system of Figure 3.3 there is a guaranteed gain margin, GM, 
given by 

GM — (3.12) 

1 ± a 

and also a guaranteed phase margin, PM, given by 


PM 


± 


-1 

cos 


[1 



(3.13) 


where the word "simultaneously" means that the gains or the phases 
of all the feedback loops may be changed at the same time within the 
limits prescribed by Equations (3.12) and (3.13) without 
destabilizing the closed-loop system. It does not mean, however , 
that the gains and phases may be changed simultaneously. 

Mukhopadhyay and Newsom [20] extended this analysis to include 
phases and gains changing simultaneously by using the particular 
form of the L matrix as given by Equation (3.7); using this form, 
the stability condition is given by 


o(L -1 - I) = /(I --J-) 2 +|- (1 - cos (p ) 
k k n 

max n n 


(3.14) 
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for all n with k n > 0. Using Equation (3.14), the general case of 
simultaneous gain and phase changes can be examined. Equation 
(3.14) can also be presented as a universal diagram for gain and 
phase margin evaluation [20]; this diagram is shown in Figure 3.4. 
For example, if the smallest _a( I + ) for a system is 0.6, then the 

closed-loop system will tolerate simultaneous gain and phase changes 
of -1.5 dB to +5.3 dB, and -30 deg to +30 deg, respectively, in all 
input loops. In a classical sense, when either gain or phase is 
changed, the margins are -4.2 dB and +8 dB or ±35 deg, 
respectively. The same results are obtained from Equations (3.12) 
and ( 3 • 1 3 ) . 

It is important to note that condition (3.10) is a conservative 
condition and that it is possible to construct a matrix, L, which 
violates the condition yet fails to destabilize the system [20] • 

This means that the predicted stability margins will be smaller than 
the actual stability margins. This is due to the fact that some of 
the small perturbations that would theoretically destabilize the 
closed-loop system will never occur in the physical system [22] • 
Nevertheless, they are still detected by a small _a(I + G-j ) - A less 
conservative evaluation of the stability margins can be obtained by 
replacing the maximum and minimum singular values in Equation (3.10) 
with the maximum and minimum magnitudes of the eigenvalues, 
respectively [20, 27]. However, in that case, the stability of the 
system can only be guaranteed if the simultaneous gain or phase 
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Figure 3.4: Universal Diagram for Multiloop Gain and Phase Margins 

Evaluation (from Reference [20]) 
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changes in all the loops are uniform [27]? i.e., all k n and <f> n are 
equal. A more detailed discussion of the use of these multiloop 
stability margins is reserved for the following chapter. 

It is important to note that the return difference matrix, 

(I + G-j ) , will have to be evaluated at each frequency of interest; 
and therefore, it is possible to construct a plot of minimum 
singular values as a function of frequency. Such a plot is called a 
sigma plot. In the next section, the characteristics of the 
singular-value analysis will be explored through the use of 
examples. 


3.2 ASSESSMENT OF THE CHARACTERISTICS OF THE 

SINGULAR VALUE ANALYSIS 

In this section, some of the basic characteristics of the 
singular-value analysis are explored and compared with the 
traditional Bode plot method. This is achieved through several 
examples that include simple third-order systems, a single- 
input/single-output (SISO) yaw damper system for a business jet, a 
multiloop yaw/roll damper, and a sampled-data SISO yaw damper for 
the business jet in Appendix B. 
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3.2.1 Simple Third-Order Systems 


It was mentioned in Section 3.1 that the singular values are a 
function of frequency and that it is therefore possible to construct 
a plot of the minimum singular value for a range of frequencies. 

Such a plot is called a sigma plot (or a-plot) and will yield a 
minimum at some frequency, a). 

Two simple, third-order systems are used to demonstrate the 
concept of the sigma plot and some of its characteristics, namely: 
how the sigma plot changes as system gain changes, how poles at high 
and low frequency affect the sigma plot, and how driving a system 
from stable to unstable affects the sigma plot. Also, the 
conservativeness of the stability margins derived from the singular 
value analysis is demonstrated. 

A) Case X. The first example system is shown in Figure 3.5. 
The system has a real pole at -2 and a pair of complex poles at 
-2 ± 4j . The frequency response is shown in Figure 3.6 for gains 
K = 40, 100, 120, and 140? while the sigma plots are shown in 
Figures 3.7 (a)-(d), at each of the different gains. The complex 
poles are unstable for K = 140. Two things are apparent. The 
minimum singular value usually does not occur at exactly the same 
frequency as the gain or phase crossover frequencies but at a 
slightly different frequency, and the frequency varies as the 
complex poles approach the imaginary axis in the s-plane. This can 
be explained by the fact that the sigma-plot analysis takes into 
account possible variations in phase and gain that are not accounted 
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Figure 3.5: Block Diagram and Root Locus of Ttiird-Order System in 

Case I 
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Figure 3.6: Frequency Response for the Third-Order System of Case I 
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for by the Bode plot. Also, note that as the gain increases and the 
stability margin of the system is reduced, the minimum singular value 
plot at the "high-frequency” range shifts down, decreasing the value 
of the minimum singular value. Therefore, since the minimum singular 
value is a measure of the closeness to the stability boundary, the 
fact that the plot shifts down means that the relative stability is 
reduced. An important observation must be made here. Note that for 
K = 140, the sigma plot has "bounced" up again. This is so because 
the minimum singular value is a measure not of absolute stability of 
a system but of how far the system is from crossing the j co-axis in 
the s-plane; i.e., how close it is to neutral stability. Also, by 
definition, ^o(A) >_ 0 (Section 3.1). 

In the "low-frequency" range, the sigma plot is shifted upward 
as K increases; this shift is related to the real pole moving to the 
left (away from the j co-axis) in the s-plane as K increases. 

The phase and gain margins from the frequency response and the 
predicted phase and gain margins from the singular value analysis 
are shown in Table 3.1. The gain margin is defined as the smallest 
gain change that will destabilize the system when the phase is kept 
constant; the phase margin is 'defined as the smallest phase change 
that will destabilize the system when the gain is kept constant. 

All the phase margins (PM*s) predicted by the singular values are 
conservative as compared to the PM from the Bode plot. The gain 
margins (GM's) from the sigma plots are also conservative except for 
gain K = 140. However, in this case, the fact that the GM from the 
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Table 3.1s Comparison of Stability Margins 
from the Bode and Sigma-Plot Analyses, 
for Case I 



Feedback 

Gain 

K = 40 

K = 100 

K = 120 

K = 140 


From Bode Analysis 
GM, dB 

to 10.35 

to 2.38 

-® to 0.804 

-0.54 to 00 


PM, deg 

00 

15.37° 

4.46° 

-7.0° 

to 

From a-Plot Analysis 





o\ 

_o 

0.61 1 2 

0.159 

0.057 

0.070 


GM, dB 

-4.15 to 8.19 

-1 .28 to 1 .508 

-0.477 to 0.505 

-0.58 to 0.649 


PM, deg 

±35.6 

±9.14° 

±3.23° 

±4.0° 



sigma-plot analysis is larger than the GM from the Bode plot is due 
to lack of resolution in the computation of the singular values; 
i.e., if the frequency interval used for the calculation is reduced, 
the situation would be reversed. In any case, the GM's calculated 
from the sigma-plot analysis are not very conservative (i.e., they 
are close to the GM's from the Bode plot) for this particular 
example. 

B) Case II. The second example is shown in Figure 3,8. This 
example is similar to Case I except that the real pole is unstable 
at +2 rather than stable at -2. The frequency response is shown in 
Figure 3,9 for gains K = 37, 50, 60, and 70. Note that the phase is 
equal to -180 deg at a frequency of zero and then again at a 
frequency between 3 and 4 rads/sec. This means that there are two 
GM's that we need to look at: one at each of the two phase- 

crossover frequencies. At K = 37, the real pole is unstable while 
the complex poles are stable; at K = 50, all the poles are stable; 
at K = 70, the real pole has become stable while the complex poles 
are now unstable. 

The sigma plots are shown in Figures 3.10(a)-(d). At K = 37, 
the minimum singular value is provided by the " low-frequency" pole 
(Figure 3.10(a)), while at K = 50 the minimum singular value is 
provided by the "high-frequency" pole (Figure 3.10(b)). For K = 60, 
the "low-frequency" range of the curve continues its shift upward, 
while the "high-frequency" part has dipped down even more, meaning 
that the distance of the real pole from the jw-axis on the s-plane 
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Figure ^3.8: 


Block Diagram and Root Locus of Third— Order System in 
Case II 
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Figure 3.10: 


Sigma-Plots of the Third-Order System in Case n, for 
Several Feedback Gains 
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has increased while that of the complex poles has decreased. This 
trend continues for K = 70, where the '’high-frequency” part of the 
curve has bounced up after reaching a = 0 at the crossing of the 
jw-axis by the complex poles. 

Note that, if the analyst looks at the sigma plot for K = 70, 
he may erroneously conclude that the GM of the system is given by 
the minimum singular value at the "high-frequency" range of the 
curve in the direction of increasing gains. It is necessary to look 
at the closed-loop eigenvalues of the system to determine where the 
system stands from the absolute stability point of view. 

The PM and GM from the Bode analysis are compared, in Table 
3.2, with the stability margins from the sigma-plot analysis. The 
conservativeness of the sigma-plot analysis is again apparent from 
the comparison. 

From the two examples presented, the following conclusions can 
be drawn: 

1 ) The PM and GM obtained from singular values are 
conservative when compared to the Bode-plot stability margins. 

2) In simple systems as the ones presented, the effect of the 
poles at different frequencies can be detected from the sigma plots. 

3) It is not possible to determine if a system is stable by 
looking only at the sigma plots. The stability margins merely tell 
the designer how far the system (or a particular pole) is from the 
stability boundary. 
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Table 3.2: Comparison of Stability Margins 

from the Bode and Sigma-Plot Analyses, 
for Case II 


Feedback 

Gain K = 37 K = 50 K = 60 


K = 70 


From Bode Analysis 





GM, dB 

-« to 0.677 

-» to 2.16 

— oo to 0.57 

-0.765 to 

PM, deg 

OO 

18.25° 

6.50° 

-7.63° 

From a- Plot Analysis 





_a 

0.075 

0.196 

0.0547 

0.0802 

GM, dB 

-0.63 to 0.67 

-1 .55 to 1 .89 

-0.463 to 0.489 

-0.67 to i 

PM, deg 

±4.29° 

±1 1 .24° 

±4.59° 




4) The sigma plot information must be complemented by the 
determination of the closed-loop eigenvalues. If the system is 
SISO, the root locus is of great help. 


3.2.2 Single-Input/Single-Output Yaw Damper for a Typical 
Business Jet 

Yaw damping stability augmentation system (SAS) is employed 

when the dutch roll mode of an airplane is poorly damped and when an 

analysis of the stability derivatives shows that increasing C is 

n r 

the most effective way to improve dutch-roll damping [31]. This is 
the case with many business and general aviation airplanes, 
particularly on final approach. 

To prevent the yaw damper from interfering with the pilot in 
turn entries and in steady turns, a washout circuit is usually added 
[31]. The washout filter assures that the rate-gyro output is 
effective as a feedback signal only during transient motions, by 
filtering out low-frequency signals and allowing the passage of only 
high frequency signals. 

The yaw damper system presented in this section was obtained 
from Reference [31]. The airplane aerodynamic and inertial 
characteristics, and the flight condition data needed for the yaw 
damper design, are those of the business jet contained in Appendix B 
of this report. The handling qualities requirements for a typical 
business jet on dutch-roll damping, and time to double amplitude for • 
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the spiral mode, are shown in Table 3.3. The open-loop airframe 
characteristics are also shown in Table 3.3, along with the open- 
loop eigenvalues. Note that the airframe dutch-roll damping is very 
low; also, although the spiral mode is slightly unstable, this mode 
is of no concern, due to the large time to double amplitude. 

The root locus and block diagram for the yaw damper are shown 

in Figure 3.11. A gain of K = 1.3 is selected to maximize the 

dutch-roll damping. The system has a positive feedback, due to the 

negative gain in the r/S transfer function; and therefore the ±0 

R 

deg criterion applies when determining the GM from the frequency 
response plot shown in Figure 3.12. The system is unstable in the 
spiral mode; but without looking at the root locus or closed-loop 
eigenvalues, it is not possible to determine this from the Bode 
plot. The GM and PM for the dutch-roll mode are 
GM = -°° to 10.6 dB 
PM = -97° and 50.4°. 

The sigma plot for this system is shown in Figure 3.13, which 
gives GM and PM of 

GM = -4.3 to 9.03 dB 
PM = ±38° 

for o_ = 0.6466. The GM obtained from the sigma plot analysis is 
very conservative in the negative direction (decreasing gain) and 
less conservative in the positive (increasing gain) direction. The 
PM appears to be very conservative in both directions. Note that 
for low frequencies (u) < 0.02 rad/sec) the sigma plot approaches 
o_ = 1 asymptotically. This means that for the low-frequency pole 
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Table 3.3: Handling Quality Requirements and Open-Loop Airframe Characteristics 

for a Typical Business Jet 



Flight Condition: 

Final Approach; 

Cat. C, Class II, 

Level 1 



Minimum 

Time to Double 
Amplitude, T~ 

Minimum 
Dutch Roll 
Damping, C D 

Minimum 

Minimum 

“b 

Maximum 
. Roll Time 
Constant T R 

Handling Quality 
Requirement 

20 secs 

0.08 

0.15 rad/sec 

0.4 rad/sec 

1 .4 sec 

Open-Loop 

Characteristics 

185.2 secs 

0.00239 

0.00263 rad/sec 

1.10 rad/sec 

1 .20 sec 


Open-Loop Eigenvalues 


Dutch-roll poles: -0.00263 ± 1.1 01 j 

Roll pole: -0.83178 

Spiral pole: +0.003743 




Yaw- Rate-To- Rudder Transfer Function: 



= - 0.708 


(s+0.912) 


•136) 2 + (.521) 2 ] 


(s + 0.832) (s - .004) [(s + .0031)2 + (1.10)2] 


Figure 3.11: Yaw Damper Root Locus for a Business Jet (taken from 

Reference [31 ] ) 
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3.13: Singular-Value Elot for Yaw Damper with Washout Filter 

(.4®/ ( 4s + 1 ) J for a Typical Business Jet, k = 1.3 
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the GM is about -6 dB and 00 in the decreasing and increasing gain 


directions, respectively, as determined from Figure 3.4. 

A second yaw-damper design was tested in which the washout 
filter was replaced by a lead-lag filter to allow the spiral root to 
move into the stable region of the s-plane. In this manner, it was 
possible to compare an unstable system with a stable one. The lead- 
lag filter used was (4s + 0.04)/(4s + 1); the Bode plot for gain K = 
1.3 is shown in Figure 3.14. Observe that, in this case, the low- 
frequency part of the curve crosses over the 0 dB line, indicating 
that the spiral root has crossed to the left-hand side of the s~ 
plane for the particular gain chosen (K = 1.3). 

At a high frequency (dutch-role poles) there is not much 
difference between the design with the lead-lag filter and the one 
with the washout filter. The GM and PM for the dutch-roll mode are 
GM = - to to 10.2 dB 
PM « 46.8° to -99°. 

For the spiral mode, the margins are 
GM = -8.4 dB to «> 

PM = 212. 

The sigma plot (Figure 3.15) yields a o_ = 0.6523 for the high- 
frequency poles, which corresponds to 
GM = -4.36 to 9.18 dB 
PM = ±38.07° 

and a "relative minimum" a = 1.587 at low frequency 

— R 

( to » 0.0 rad/sec), which gives 
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Figure 3.14: Frequency Response for a Yaw Damper with Lead-Lag 

Filter (4s + 0.04) /(4s + 1) for a Typical Business 
Jet, K = 1 .3 
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Figure 3.15: 


Singular-Value plot for a Yaw-Damper with 
Filter (4s + 0.04)/(4s + 1 ) f or a Typical 
Jet, K = 1.3 


Lead-Lag 

Business 
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GM = -8 • 25 dB to « 


PM = ±105°. 

In both cases the stability margins obtained from the sigma 
plot analysis are conservative with respect to the Bode plot 
analysis. Note that the low-frequency part of the sigma plot has 
shifted upward in relation to the plot in Figure 3.13; it is evident 
that the sigma plot alone is not enough to understand what is 
happening with the system. 

To demonstrate the effects of the spiral root more clearly, the 

gain was changed from K = 1.3 to K = 0.65 for the yaw-damper design 

with lead-lag filter. The gain change reduced the GM for the low- 

frequency range of the curve, as shown in the Bode plot of Figure 

3.16. The same effect is detected in the sigma plot (Figure 3.17) 

by a downward movement of the low-frequency range of the curve. In 

this case, the minimum _a has become the one at low frequency with 

£ = 0.2932, while the former minimum has become a "relative-minimum" 

with a « 0.794 • Therefore, at K = 0.65, the spiral pole is closer 
R 

to instability than the dutch-roll poles. 

In the following section, the yaw-damper design is augmented 
with a roll damper to improve the dutch- roll damping beyond what was 
achieved with the yaw-damper system. This arrangement also provides 
the opportunity to analyze a multi loop control system. 
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Figure 3.17: Singular-Value Plot for a Yaw-Damper with Lead-Lag 

Filter (4s + 0.04)/(4s + 1) for a Typical Business 
Jet, K = 0.65 
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3*2*3 A Multi loop Yaw/Roll Damper System for a Typical Business Jet 


The yaw/roll-damper system for the business jet in Appendix B 

is shown in Figure 3.18. A plot showing how the locus for the 

dutch-roll poles changes for several values of roll-damper-loop 

gain, K* , is contained in Figure 3.19. A gain of Kx =1.0 for the 
A °R 

yaw-damper loop, and a gain of Kr =0.2 for the roll-damper loop 

°A 

were chosen for the multi loop design. 

Traditionally, the multiloop design would be analyzed using a 
frequency response method (example: Bode plot) following a one- 

loop-at-a-time approach. This situation is depicted schematically 
in Figure 3.20. The GM and PM of each loop are each obtained 
independently, by "closing" the other loop with its respective 
design values and reducing it as part of the dynamics of the loop 
being analyzed. In this manner, the stability margins of either 
loop can be estimated when that particular loop is subjected to 
uncertainties, and assuming that the parameters of the other loop 
are kept constant. The problem with this approach is that the 
analysis will not give any information on the stability margins of 
the whole system when subjected to simultaneous gain or phase 
uncertainties in both loops. It can be expected that the stability 
margins for each loop are dependent, to some degree, on the actual 
parameters of the other loop. This is a case in which the use of 
the sigma-plot analysis is useful in determining the multiloop 
stability margins. 
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In order to explore how the sigma-plot analysis compares to the 
traditional one-loop-at-a-time Bode analysis, the following analyses 
were made : 

1 ) Frequency response of the yaw-damper loop with the roll- 
damper loop "closed" at the design gain value K* =0.2. 

2) Sigma-plot analysis for case (1) above. Note that this is 
not the multiloop analysis but a one-loop-at-a-time sigma-plot 
analysis used to explore how the Bode and sigma-plot analyses 
compare . 

3) Frequency response of the roll-damper loop with the yaw- 

damper loop closed at the design gain value = 1.0. 

R 

4) Sigma-plot anlaysis for case (3) above. 

5) Multiloop sigma-plot analysis. 

6) Eigenvalue-plot analysis. 

In the eigenvalue-plot analysis, the minimum magnitude of the 
eigenvalues of the return difference matrix, (I + G 1 ) , is plotted 
along with the minimum singular values for a range of frequencies. 

By substituting the minimum and maximum magnitudes of the 
eigenvalues for the minimum and maximum singular values respectively 
in Equation (3.10), less conservative PM and GM can be obtained 
[27]. 

The frequency response analysis for the yaw-damper loop in 
Figure 3.21 (±0° criterion) gives 

GM = -°° to 21 dB at to = 0.49 rad/sec 
PM = +90° at a) = 0.87 rad/sec 

-108.4° at 03 = 1.42 rad/sec. 
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Figure 3*21: Frequency Response of the Yaw-Damper Loop of a 

Yaw/Roll Damper Multiloop System for =1.0. 

R 

(This "single-loop" analysis includes the roll-damper 

loop, closed at Kr = 0.2, in the dynamics of the 
°A 

yaw-damper loop. ) 


The corresponding sigma plot in Figure 3.22 shows 


o_ = 0.8862 at 0 ) = 0.547 rad/sec 
GM = -5.51 to 18.87 dB 
PM = ±52.6°. 

The stability margins from the sigma-plot analysis are 
conservative compared to the frequency response analysis: the GM is 

very conservative in the direction of decreasing gains, and the PM 
is very conservative in both directions. Note that the frequency 
for minimum o_ does not match any of the crossover frequencies but 
lies at a slightly different frequency; this is usually the case 
[20, 34] and can be explained by the fact that the singular-value 
analysis accounts for simultaneous changes in phase and gain. If we 
look at the frequency response plot, at the frequency for minimum 
a ( oj = 0.547) the worst direction of simultaneous gain and phase 
change is found to be about 18 dB and 38° [20, 34] . If now we look 
at the phase crossover frequency (a) = 0.49 rad/sec) on the sigma 
plot, we find that the predicted GM is 21 dB, matching the GM from 
the Bode analysis [20]. 

Similarly for the roll-damper loop, the frequency reponse in 
Figure 3.23 (±180° criterion) shows 
GM = ±» 

PM = ±°° 

while the sigma-plot analysis in Figure 3.24 shows 
o_ = 0.97237 

GM = -5.89 to 31 .17 dB 
PM = ±58°. 
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MINIMUM SINGULAR VALUE, 




Loop 


FREQUENCY, o) (rad/ sec) 


sncy Response for Roll-Damper Loop of a Yaw/Roll 
: Multi loop System, for K* =1.0 and Roll-Damper 

°R 

5ain, K 5 = o.2 


Similar conclusions can be drawn for this analysis as for the 


yaw-damper loop analysis , and they will not be repeated here. 

The next analysis of interest is given by the multiloop 
singular value plot. Figure 3.25 shows the sigma plot (lower curve) 
and the eigenvalue plot. Two things are apparent from the multi loop 
sigma plot: the minimum o_ is much smaller than either of the two 

previous single-loop minimum o_ of Figures 3.22 and 3.24. This means 
that, according to the multiloop analysis, the stability margins are 
smaller than the ones obtained through the one- loop-at-a- time 
approach. The multi loop margins obtained are 
GM = -3.53 to 6.05 dB 
PM = ±29.1 ° 

for o_ = 0.502 at go = 0.76 rad/sec. 

The important characteristic of these multiloop margins is that . 

they represent the limits within which the gains or the phases of 

the two feedback loops may be changed simultaneously and 

independently without destabilizing the closed-loop system. In 

other words, Kr can be changed by +5 dB while Kr is changed by 
k °A 

-3 dB, and the system would still be stable. Furthermore, if we 
were interested to know the permissible range for gain variations 
considering that the phases may simultaneously change within ±20°, 
Figure 3.4 could be used to obtain a GM = -2.2 to 4.8 dB. 

The upper curve in Figure 3.25 is a frequency plot of the 
minimum magnitude of the eigenvalues of (I + G^), which represents 
an upper boundary of the _a(l + ) plot [20, 21]. The 
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FREQUENCY, to (rad /sec) 


Figure 3*24: Singular-Value Plot for Roll-Damper Loop of a Yaw/Roll 

Damper Multiloop System, for Kx =0,2 and Yaw-Damper 

°A 

Loop, Kx =1.0 
°R 



Figure 3*25: Multiloop Singular— Value and Eigenvalue Plot for 

Yaw/Roll Damper System with Yaw-Damper Loop Gain, 

K 6 = 1#0 ' and Roll-Damper Loop Kx = 0.2 
R Ojv 
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min|A^(I + ) | represents the radial distance of the eigenvalue 

closest to the origin and is a measure of the closeness of (I + 
to singularity in a limited sense; i.e., for uniform perturbations 
of gain or phase in each of the loops. 

By uniform perturbations, it should be understood that if the 
perturbation matrix, L, is given as in Equation (3.7) 


L = 

and if the gain change in one loop is k j , then the gain changes in 
the remaining loops must also be given by k^ = kj , for i = 1,n. The 
same restriction applies to the simultaneous variations in phase in 
each loop [27]. Yeh, et al. [27], claim that this is the least 
conservative computation that is possible when a norm-bounded robust 
stability criterion is used. 

Mukhopadhyay and Newsom [20] propose the use of the eigenvalue 
plot to obtain a set of "relaxed-stability” margins, where the gain 
or phase in each loop can be allowed to change in any manner within 
the prescribed limits. The only shortcoming of these limits is that 
the stability of the closed-loop system is not -always guaranteed, 
unless the phase or gain variations are uniform in the sense 
described by Yeh, et al. [27]. 




V 


k 2 6 


j 4>- 


(3.7) 


j <t> 


k e 
n 
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Independently of how the limits obtained from the eigenvalue 
plot are interpreted (uniform or nonuniform), they are given by 
GM = -5.29 to 15.88 db 


PM = ±49.6° 

for min | | = 0.839 at w = 0.566 rad/sec. 

Observe that these margins are considerably less conservative 

than the ones obtained from the multiloop sigma plot, although still 

conservative compared to the Bode analysis. 

The plot in Figure 3.26 contains the actual stability boundary 

for the dutch-roll mode of the airplane, for combinations of gains 

in the yaw- and roll-damper loops, and assuming no variations in 

phase. For a simple example like the one presented here, this 

stability boundary can be calculated very easily by an iterative 

procedure. For a more complex problem, the calculation of the 

actual stability boundary is very cumbersome; however, this example 

allows us to show a very interesting and important situation. 

The design point is given by Kr =1.0 and K* = 0.2. The 

°R °A 

stability boundary given by the ABCD block contains the combinatins 

of Kg and Kg which, according to the Bode analysis, will not 

destabilize the system (point D is not shown in the figure; it lies 

at Kf = 00 and K* =0.5. This situation is also present in the 
°A °R 

other boundary blocks to avoid overcrowding the figure.) Note that. 


for the large amount of possible combinations of Kf and 


Kg that 


lie to the right of the curve, the dutch-roll mode will be 


unstable. This situation is not detected by doing a one-loop-at-a- 


time Bode analysis. 
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FEEDBACK CAIN FOR THE ROLL DAMPER LOOP 



Figure 3.26: Comparison of Predicted and Actual Stability 

boundaries for the Dutch-Roll Poles of a Business Jet. 
(Block MNOP represents the boundary predicted by the 
multiloop singular-value analysis.) 
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Boundary EFGH (point H is not shown) is given by the "SISO 


singular-value analysis. The amount of unstable combinations of 

R 

and Kr that have not been predicted has been reduced considerably 
°A 

due to the conservativeness of the margins; however, since it does 

not account for simultaneous variations of gains in both loops, 

stability is not always guaranteed. 

Boundary IJKL (point L is not shown) is given by the "relaxed- 

stability" margins obtained using the eigenvalue plot instead of the 

sigma plot. Although stability is generally not guaranteed unless 

K* and Kf change by the same factor, for our present example there 
°R °A 

are no possible unstable combinations of and Kr . Therefore, in 

°R °A 

this case, the relaxed-stability margins seem to provide the least 
conservative approach to the determination of the multiloop 
stability margins. However, it is very important to remember that 
we have assumed that there is no change in the phase of each loop. 

If we considered simultaneous phase and gain changes, the limits 


would be reduced. 


The most conservative of the limits are given by the multiloop 
sigma-plot analysis and are represented by block MNOP (point P is 
not shown). Although very conservative, these limits guarantee that 
the closed loop system will be stable when the gains are varied, in 
any manner, within the limits. In this case it has also been 
assumed that there is no change in phase; otherwise, the limits 
would be reduced further. 
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3.2.4 A Sampled-Data Yaw-Damper System for a Typical Business Jet 


This section explores briefly the application of the singular- 

value concepts to a sampled-data yaw-damper SAS. The system is 

represented in the block diagram of Figure 3.27. The transfer 

function, r(s)/6 (s) , is the same as for the continuous system 
R 

given in Figure 3.11. 

The fundamental concepts of the singular value analysis 
presented in Section 3.1 are also applicable to sampled data 
systems; i.e., the stability condition is given by 

£(I + G ( z ) ) > 0 (3.15) 

where G^z) is evaluated at 

jwT AT \ 

z = e J (3.16) 

where z is the z- transform variable, and T is the sampling period 

[25]. However, the concepts of simultaneous phase and gain margins 

(Equations 3.12 and 3.13) were derived for continuous systems and 

there is no information in the literature relative to their 

applicability to sampled-data and discrete systems. 

In the following analysis, no attempt was made to prove (or 

disprove) theoretically the applicability of Equations (3.12) and 

(3.13) to sampled-data systems. Instead, the concept of stability 

margins was explored by predicting the phase and gain margins, of a 

sampled-data yaw-damper SAS, using singular values and then 

comparing them with the phase and gain margins obtained from a 

frequency response in the w* -plane. 
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The transformation from the z-plane to the w' -plane is made by 
using [38] 


, 2 z - 1 

w» = — 


(3.17) 


T z + 1 

This transformation is used to regain the resolution lost as z takes 
on values around the unit circle. Also the use of the w* -plane 
allows us to regain the simplicity of the Bode frequency response 
technique which is not present in the z-domain due to the 
nonrational nature of the discrete transfer functions, where the 
frequency appears in the form z = e^ a)T . 

The w* -transform of Equation (3.17) can also be expressed as 


w' = — 


2 e sT - 1 2 


sT 


T sT 
e + 1 


= — tan h — 
T 2 


(3.18) 


and, for the frequency response (s = j 03 ) we have 


. 2 


ojT 


w* = jv = j - tan — 


(3.19) 


T "2 

and we see that while z goes around the unit circle, the w* -plane 

frequency, v, stays real and goes from 0 to 

The frequency response for the sampled-data yaw-damper SAS is 

shown in Figure 3.28 for a design gain of K* =0.39 and a sampling 

°R 

period of T = 0.8 seconds. The corresponding singular value plot is 
shown in Figure 3.29. 

Note that the Bode plot is calculated in the w* -plane, and 
therefore the frequency in the horizontal axis of Figure 3.28 is the 
w* -plane frequency, v. Although the singular value analysis is done 
using the s-plane frequency, a), the horizontal axis in Figure 3.29 
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W' -PLANE FREQUENCY, V (rad/sec) 


Figure 3.28: Frequency Response for Sampled-Data Yaw-Damper SAS 

with K* 

6 R - 


0.39 and T = 0.3 Seconds 



W' -PLANE FREQUENCY, V (rad/sec) 


Figure 3.29: Singular-Value Plot for Sarapled-Data 

with K* _ 0.39 and T = 0.3 Seconds 


Yaw-Damper 


also shows the w 1 -plane frequency, v. This transformation was done 
using Equation (3.19) with the objective of having the Bode and 
sigma plots expressed as functions of the same frequency. 

The stability margins from the frequency response computation 
are given by (±0° criterion) 

-«> to 9.53 dB @ v = 0.5575 rad/sec 
(or 0 ) = 0.5562 rad/sec) 


G M = 


-» to 11.52 dB 0 v = 4.107 rad/sec 
(or to = 3.68 rad/sec) 

62.72° 0 v = 0.8319 rad/sec 
(or to = 0.8277 rad/sec) 

PM = 

-68.4° 0 v = 1.6164 rad/sec 
(or to = 1 .586 rad/sec) 

The predicted stability margins from the sigma plot (Figure 
3.29) are 


GM = -4.35 to 9.12 dB 
PM = ±37.9° 

for a_ = 0.65007 § to = 0.6015 rad/sec 
(v = 0.6032 rad/sec) 

As it occurrs in the continuous-system example (Section 3.2.2), 
the GM predicted from the singular value computation is very 
conservative in the direction of decreasing gains and only slightly 
conservative in the direction of increasing gains. The PM predicted 
is just as conservative as predicted for the continuous system in 
Section 3.2.2. In general, the relative stability of this sampled- 
data design is very similar to the relative stability of the 
continuous-system design of Section 3.2.2. 
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An important check on the correspondence between the stability 
margins from the Bode analysis in the w 1 -plane, and the stability 
margins predicted from the singular value analysis, is given by the 
fact that [20], for a SISO system, the stability margins predicted 
by the minimum singular values at the crossover frequencies (phase- 
and-gain-crossover frequencies) must match the stability margins 
from the Bode analysis. This correspondence is shown below to exist 
for the sampled-data yaw-damper SAS analyzed here: 

Phase-Crossover GM from Predicted GM from 

Frequency, v Bode Plot Sigma Plot 

(rad/sec) (dB) (dB) 

0.5575 to 9.53 -4.46 to 9.65 

4.107 -* to 11.52 -4.78 to 11.54 

The match of GM in the increasing-gain direction is achieved 

within the accuracy permitted by the analysis. Note that no match 

exists for the GM in the decreasing-gain direction. This is also 

true for continuous systems, and we do not have an explanation for 

this occurrence at the present time. In the case of the PM's, the 

analysis shows 

Gain-Crossover PM from Predicted PM from 

Frequency, v Bode Plot Sigma Plot 

(rad/sec) (degrees) (degrees) 

0.9319 62.72 ±63.05 

1.6164 -68.4 ±68.4 

Again, the match exists at least in one direction (decreasing or 

increasing phase). 
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From the results presented above, we feel that there is 
correspondence between the stability margins obtained from the w'- 
plane frequency response and those predicted from the sigma plot, 
for sampled-data systems. However, we suggest that further research 
should be carried out to prove the validity of Equations (3.12) and 
(3.13) when used in analyzing sampled-data systems. Also the 
conditions for such validity should be investigated. 

Some important conclusions about the use of the multiloop 
singular-value analysis are as follows: 

1 ) The use of the one-loop-at-a-time frequency response 
methods to analyze multiloop control systems does not allow the 
analyst always to detect instabilities which arise when simultaneous 
changes in gain or phase in each loop are considered. 

2) The stability margins obtained from the singular-value 
analysis are very conservative in some directions compared to the 
limits obtained through the Bode-plot analysis; this applies to 
multiloop and SISO systems. 

3) The use of the limits obtained from both the sigma- and 
eigenvalue-plot analysis seems to provide the best alternative, 
today, to the use of the SISO Bode plot in the analysis of multiloop 
control systems. The sigma plot provides a lower boundary while the 
eigenvalue provides the upper boundary to the relative stability of 
the multiloop systems. If these two plots are close together, the 
conservativeness of the analysis is small. 
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4) The use of the singular values in a way similar to the 
one-loop-at-a-time frequency response, in the analysis of multiloop 
systems, also provides some insight for potential problems that 
could arise with simultaneous changes in phase and gain in each 
loop. 

5) The use of singular-value analysis to obtain stability 
margins of sampled-data systems seems feasible; however, more 
research is recommended in this area. 

In this chapter we have presented the fundamental concepts of 
singular values and several SISO and MIMO control system examples to 
explore the use of the singular-value analysis and how it compares 
with the traditional Bode method. In the following chapter we 
present the derivation of the singular-value gradients and discuss 
the use of the sensitivity analysis technique developed. A simple 
third-order example is used to show how to interpret the informaton 
provided by the sensitivity analysis, and the multiloop yaw/roll 
damper system for the business jet is used to demonstrate how the 
method can be used in determining the important model parameters. 
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4. MULTILOOP SENSITIVITY ANALYSIS METHOD 


In this chapter, expressions for the singular-value gradients 
are developed, and their use in a sensitivity analysis technique is 
discussed. The application of this sensitivity analysis technique 
is demonstrated with a simple third-order example and the multi loop 
yaw/roll damper system of the business jet in Appendix B. 


4.1 MATHEMATICAL DEVELOPMENT OF THE SINGULAR-VALUE GRADIENTS 

In order to determine the effect that changes in the system and 
controller parameters have on the relative stability of a closed- 
loop system, the derivatives of the singular values with respect to 
those parameters must be evaluated. Newsom and Mukhopadhyay [24] 
determined the gradients of the singular value_a(I + ) , with 

respect to controller parameters, and used them to perform a 
numerical optimization to search for the controller design variables 
that increase the minimum singular value of the system return 
difference matrix. 

For the present application what is desired is a measure of the 
effect on system relative stability caused by variations in the 
system and controller models. Therefore, the gradients of the 
singular value with respect to elements in the system and controller 
matrices are required. 
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Consider the state-space representation of the control system 
of Figure 3.1 to be given by 
Plant: 

X = AX + BO (4.1) 

Z = TX (4.2) 

Control Law: 

U = -KZ + R (4.3) 

Equation (4.1) represents a plant of order Ns having No output 
measurements, Z, modeled by Equation (4.2) and Nc control inputs, 

U. Equation (4.3) represents the feedback control law driven by the 
sensor output, Z, and reference input signal, R. In terms of 
transfer matrices (taking Laplace transforms), the control law is 
given by 


U(s) = -K[T(Is - A) -1 B]D(s) + R(s) 
Therefore, the control input can be written as 

u(s) = (i + H(s)G(s) ) -1 R(s) 
where K = H(s) and G(s) = [T(Is - A) -1 B]. 


(4.4) 


(4.5) 


Equations (4.1) to (4.3) can be written in an augmented form as 


X = AX + BU 


(4.6) 


U = -CX + R 


(4.7) 


where a = A, 5 = B, and C = KT. Then the Laplace transform of the 
control input can be expressed as follows (the argument " (s) M has 
been omitted) : 


O = [I + C( Is - A) -1 B] -1 R 


(4.8) 
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and, therefore, the return difference matrtix, (I + G 1 ) , can also be 
represented as 

(I + G 1 ) = (i + C( Is - A) _1 b] (4.9) 

In the case in which the control law includes controller 
dynamics, a similar derivation would be involved in the computation 
of A, B, and C matrices augmented with the controller dynamics. 

The singular values of (I + ) are a^, and the corresponding 

right and left normalized eigenvectors are and u^, respectively 
[20, 24]. Hence by definition 

(I + )v i = u i a i (4.10) 

(I + G ) u. = v.a. (4.1 1 ) 

1 i ii 

* 

for i = 1, 2, ... Nc, where A means the conjugate transpose of A. 

The normalized eigenvectors satisfy the following orthogonal 
properties : 

u. u. = 6. . and v. v. = 6. . (4.12) 

1 3 1 1 

where 6— is the Kroneker delta which is unity when i = j and zero 
when i ^ j . 

Let p be a parameter for which sensitivity information is 
needed. Differentiating Equations (4.10) and (4.11) with respect to 

if ^ 

p and then premultiplying the result by and v^ , respectively 
(where a* is the conjugate transpose of a), and adding them 
together, one obtains 
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* 3(1 + G ) . 3(1 + G ) . . 9v 

U. r V. + V. U. +(u. (I + G ) - V. a.) -r + 

X 3p XI dp XX 1 X X dp 

^ ^ u ■ 8a. ^ ^ 

+ (v. (I + G ) - u. a. ) *r — ~ = -r— (u. U. + V. V. ) (4.13) 

i 1 x i 8p 8p x x x x 


Using Equations (4.10) to (4.12) in (4.13), one obtains 

* 3(1 + G,)* 


30i 1 *3(I+'G 1 ) 

3T = 2 (U i “15 V i + V i 


3p 


u. ) 

X 


• (4.14) 


where the first and second terms in the right-hand side of the 
equation are complex conjugates. Therefore, Equation (4.14) can be 


written as 
3o 


L * 3(1 + G 1 ) 

- — = Real part of [u. r v. ] 

op X dp X 


= Re[u 


* 3(1 + G 1 ) 


3p 


V i ] 


(4.15) 


Using Equation (4.9) and letting 


$ = (Is - A) 1 


(4.16) 


Then, Equation (4.15) can be written as 


3o jL (I + G 1 ) , 3 ( I + C$B) *, 

35 — - = Re * tr [ — 35 v i u i ] 


(4.17) 


where tr[W] means the trace of W and is equal to the sum of the 
elements in the principal diagonal of W. Equation (4.17) can be 
expanded as 


3 o A 


' -j 37'' _ _ fill — ^ A _ * 

— = Re • tr [ {-r— $B + C$ "X— + C$ -r — $B} V. U. ] 

8p 8p dp dp XX 


(4.18) 


It is now possible to obtain three expressions for 8a^/8p, one 
each for the following cases: 
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1) The parameter p is an element of A; i.e., P = A 

2) The parameter p is an element of B; i.e., P = B 

3 The parameter p is an element of C; i.e., P = C. 

Then, we obtain 

^ a i - 3A - * 

- — = Re • tr[C$ " 7 T — $Bv.u. ] 

3p 3p x x 

A A 


3<?i __ gg * 

- — = Re • tr[C$ - — v.u. ] 
3p 3p x l 

B B 


(4.19) 


(4.20) 


ou i 3C - * 

-r — = Re • tr [-r — $Bv.u. ] (4.21) 

3p 3p x x 

C C 

By recalling the following matrix operation [30] , 

•|t- [tr{BAC}] = B T C T (4.22) 

3A 

and the matrix trace property [24], 

Re • tr ( A) = Re • tr(A*) (4.23) 

it is possible to extend Equations (4.19) through (4.21) to all the 

elements p , p , and p of the A, B, and C matrices, respectively. 

A B_ C 

Therefore, for a, 

3 * 

= — [Re • tr {c$A$B( v. u. )}] 

3A 3A 11 

= Re[(C$) (<£Bv.u. ) ] (4.24) 

x x 

or, using the transpose of A in order to avoid transposing the 
matrices within brackets. 
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(4.25) 



(4.26) 


(4.27) 


(4.28) 


(4.29) 


Expressions (4.25), (4.27), and (4.29) can be used to evaluate 
the singular- value gradients with respect to elements of the system 
and controller matrices. Note that the gradients, like the singular 
values, are also a function of frequency; and therefore a singular- 
value-gradients plot (or sigma-gradients plot) can be obtained for a 
range of frequencies. Also, note that the information necessary to 
obtain the gradients is evaluated in the calculation of the 
singular-values plot (or sigma plot), since $ has to be evaluated at 
each frequency, and Vj_ and are a product of the singular value 
decomposition routine used in the calculation of the singular 
values. Therefore, little additional computational effort is needed 


73 



to calculate the sigma gradients for each singular value in a range 
of frequencies, beyond what is required to evaluate the singular 
values. 

In the following section, the information from the sigma plot 
and sigma-gradients plot is combined into a sensitivity analysis 
method for multiloop control systems. 

4.2 SENSITIVITY ANALYSIS METHOD FOR MULTILOOP SYSTEMS 

There are three important aspects of a sensitivity analysis of 
a control system that need to be considered so that the information 
provided by the analysis is useful. The first is a measure of how 
sensitive the particular control system is to changes in the 
elements of the system model and controller. The second is some 
information (or engineering judgement) on the magnitude of the 
expected changes in the different system and controller parameters, 
according to the degree of refinement of the models. The third 
point that needs to be considered is the effect of those changes on 
the relative stability of the closed-loop system. In this manner, a 
control system may be very sensitive to inaccuracies in a certain 
parameter, p^; but if the model is considered to be very accurate at 
the time of the analysis, the effect of the expected inaccuracy on 
the relative stability of the system may be negligible. On the 
other hand, if the sensitivity of the system to changes in p^ is not 
very high but the expected change in p^ is large, the closed-loop 
system may be driven unstable. A reduction in the stability margins 
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may be very important for a certain closed-loop system? but if the 
nominal system was very stable, the reduction may not be very 
relevant. These examples illustrate that a sensitivity analysis 
needs to include enough information to determine what parameter 
changes are most important. The method presented here does not 
attempt to obtain a "one-number" answer but, rather, to combine 
several tools to help the designer (or analyst) determine which 
parameters are most important for his (her) particular problem. 

The sensitivity analysis method follows the steps illustrated 
in Figure 4.1. First, the sigma plot of the nominal control system 
is calculated, along with the closed-loop system eigenvalues. The 
sigma plot will tell the designer what is the minimum singular 
value, and therefore the phase and gain margins of the nominal 
system can be evaluated using the universal diagram (Figure 3.4) or 
Equations (3.12) and (3.13). If the nominal system is stable, the 
phase and gain margins will only determine "how far" the system is 
from the stability boundary. If the system is unstable, the phase 
and gain margins will again determine "how far" the system is from 
the stability boundary; however, they will not give information on 
the absolute stability of the system. (The same kind of information 
is obtained from a Bode Plot in the SISO case.) Thus, the need for 
evaluating the eigenvalues of the closed- loop system to determine 
absolute stability. It must be stressed that the closed-loop 
eigenvalues must not be confused with the eigenvalues of the return 
difference matrix, which were mentioned in Section 3.2 as an , 
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Figure 4.1: Multiloop Sensitivity Analysis Technique 
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additional measure for judging relative stability. The second step 
in the sensitivity analysis is to evaluate the singular-value 
gradients (sigma-gradients plots) for the nominal closed-loop 
system. This information will tell the designer which parameter 
variations (or inaccuracies) are the most important from the 
relative stability point of view. With the aid of information about 
the expected variations of those parameters, the third step is to 
select the parameters of interest. This also includes a look at the 
values of the sigma gradients at the different frequencies and how 
they compare with the minimum singular value at that frequency. The 
next step is to evaluate the sigma plot for the closed-loop system 
perturbed by the variations on the selected parameters. Comparing 
the nominal-system and perturbed- system sigma plots, the designer 
can judge the overall effect of those variations on the relative 
stability of the system, and check if the correct parameters were 
selected. 

An illustration of the use of the sensitivity analysis is 
included in Section 4.3, using a simple third-order system and a 
yaw/roll damper system for a typical business jet. 


4.3 APPLICATION OF THE SENSITIVITY ANALYSYS TECHNIQUE 

In this section, the sensitivity analysis technique is applied 
to two different feedback control systems. The first example is a 
SISO third-order system; the objective is to show how the technique 
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can be used and the information interpreted. The second example is 
the multiloop yaw/roll damper design presented in Section 3.23; the 
objective in this case is to show how the sensitivity analysis 
information can be used to select the model parameters that have the 
greatest effect on the relative stability of the system. As we have 
discussed before, this information can be very useful in the 
selection and design of flight and wind-tunnel tests, when the 
objective is the optimization of the allocation of such resources as 
budget, flying time, wind-tunnel time, and instrumentation. 


4.3*1 Sensitivity Analysis of a Simple Tftird-Order System 

The single-input/single-output system we are dealing with in 
this section is shown in Figure 4.2. We will carry out a 
sensitivity analysis for two cases. Case I explores the sensitivity 
of the system to an error in the location of the zero? i.e., when 
the plant transfer function is given by 

G(s) = — j (4.30) 

(s + 2) (s+ 4s + 20) 

and in Case II we look at the sensitivity of the system to an error 
in the location of the second-order poles; i.e., when the plant 
transfer function is given by 

G(s) = ~ 2 (4.31) 

(s + 2) [ (s + B> + 4 ] 
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Figure 4 * 3 : Representation of the Plant Given by 

\ _ s + ot 

5 ' s ' "3 2 Used to Find 

s + 6s + 28 s + 40 

the State-Space Form 
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A) Case I. The state-space representation of the system is 
given by 


• 


—■ 



— 


— 

X 1 


0 

1 

0 

x i 


0 

X 2 

= 

0 

0 

1 

X 2 

+ 

0 

• 








X_ 


-40 

-28 

-6 



1 

3 


— 



3 _ 




(4.32) 


and 


z = [a 


0 ] 



(4.33) 


for the plant, where x.| = where £ is a variable used as a 
mathematical artifact in the representation of the plant (see Figure 
4.3), and z = 0; and 


6 = -Kz (4.34) 

for the feedback-control law, when 0 =0. Equations (4.32) to 

com 

(4.34) are equivalent to Equations (4.1) to (4.3). The system can 
then be reduced to the form given by Equations (4.6) and (4.7), by 
using Equation (4.32) and 


6 = -[ Ka K 


— - 


- - 




1 


1 


= -c 

x^ 

2 


2 



x^ 

L 3 _ 


_ 3_ 


(4.35) 


In this case, we are interested in 9a/ 3a ? i.e., the 
sensitivity of the relative stability of the system to changes in a; 
and it is obtained through Equation (4.29). 
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The feedback gain, K, chosen for the nominal design (a = 0) is 


K = 200; the corresponding closed-loop eigenvalues are 


pole 

real part 

imaginary part 

1 

-0.17623 

0.0 

2 

-2.91188 

14.78 

3 

-2.91188 

i 

. 

00 


The singular value plot for the nominal system is shown in 
Figure 4.4, with the minimum o_ occurring at a) = 15.81 rad/sec. The 
stability margins are 

GM = -2.89 to 4.36 dB 
PM = ±22.7° 
for a_ = 0.3948. 

Since Equation (4.29) calculates the singular-value gradients 
with respect to elements in the C matrix. Figure 4.5 shows 8a/3(Ka), 
which is the gradient with respect to element C(1,1); and Figure 4.6 
presents the singular-value gradient with respect to a, 3a/ 3a, which 
was obtained by using 
3_a 3_a 

- 5 - - K v ; t (4.36) 

3a 3(Ka) 

The singular-value gradient plots show that the system is very 
sensitive with respect to a (i.e., the gradient is large) at 
frequencies below gj = 1.0 rad/sec, with a maximum sensitivity at w = 
0. The importance of this sensitivity will alsQ depend on the 
magnitude of the perturbation (i.e., the value of a) and on the 
value of a at those frequencies. 
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Figure 4.4: Sigma Plot of the Nominal System (a = 0, 6=2) 
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FREQUENCY, oj RAD/ SEC 

Figure 4.5: Gradient of the Minimum Singular Value with Respect to 
Ka f i.e. 8a/9(Ka) , as a Function of Frequency 
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Figure 4,6: Frequency Plot of the Gradient of a 



Figure 4,7: Sigma Plot of the Perturbed System (a 

the _a at a Frequency of u « 0 rad/sec 


with a 


.2) Showing 
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If we let a = -0.2, which means that the zero is on the right 
hand side of the s-plane, the estimated change in _a, at u = 0 
rad/sec, is A o_ « -1.0; this means that the stability margins at 
very low frequencies will be considerably reduced. This is shown in 
the sigma plot of Figure 4.7, where the minimum o_ now appears at 
a) = 0 rad/sec. Also, the closed-loop eigenvalues of the perturbed 
system are 

pole real part imaginary part 


1 0.0 0.0 

2 -3.0 14.79 

3 -3.0 -14.79 


which confirms the results of the sensitivity analysis. 

B) Case IX. The state-space representation of the system is 
given by 


— 


— 


- 


- — 


- 1 

• 

X 1 


0 

1 

0 


X 1 


0 

. 

X 2 

= 

0 

0 

1 


X 2 

+ 

0 

. 

X 3 


f 1 (3) 

f 2 (e> 

f 3 <e>_ 


_ X 3_ 


1 


and 


z = [ 0 



(4.38) 


for the plant, where x 1 = c, £ is defined in Figure 4.3, and z = 0; 
and 


84 



6 = -[ 0 K 0 ] x 


(4.39) 


In this case, we are interested in 3q/3f ^ , do/df^r 3a/3f^ # and 
3a/3B» Since f-j(B), f 2 (B)r and ( 3 ) are elements of the A matrix. 
Equation (4.25) must be used. 

The closed-loop eigenvalues and the sigma-plot for the nominal 
system (K = 200, a = 0, B = 2) are the same as in Case I. The 
singular-value gradients are plotted in Figure 4.8. Note that these 
gradients are important only in the frequency range between 1.0 and 
20.0 rad/sec, with 3a/3f^ being the most important. 

In order to obtain the sensitivity with respect to B, the 
contribution of each of the elements, 3a/3f , 3a/3f , and 3a,.,. has 

— 1 — 2 — / 3 f 3 

to be accounted for. This is done by using the chain rule to obtain 


3_o 3f 1 3_a 3f 2 3a 3f 3 3_a 

T3 = IF" FT" + FiT 9f 2 + FF“ 3f 3 


(4.40) 


where 9f^/33, 9f 2 /93, and 9f 3 /93 are determined from 


f 1 (0) = -2(0 + 16) and 


= -43 


(4.41 ) 


7 ol 2 

f„( 3) = -(3 + 43 + 16) and — - = -2(3 + 2) (4.42) 

z dp 


f 3 ( 3) = -2(3 +D and = -2 


(4.43) 


which, when evaluated at nominal conditions (3 - 2), give 
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-GRADIENT, 3a/3f. 



.01 0.1 1.0 10. 100. 


FREQUENCY, w RAD/ SEC 


Figure 4.8: 


The singular-value gradients with respect to f 

and are important only between 1.0 and 20 rad/s 



3a 


3a 


3a 


(4.44) 


3a 

"33 = " 8 lf~ ~ 8 HfT ~ 2 "3f^ 

The plot for the gradient with respect to 3 is shown in Figure 
4.9. At low frequency the gradient approaches zero; and it reaches 
a maximum in the negative direction at w = 3.84 rad/sec, with 
3o_/3 B = -5.7. However, the sigma plot in Figure 4.4 shows that, at 
that frequency, the o_ = 11.7; and, therefore, the effect of any 
expected changes in 3 will not be very important. Figure 4.9 also 
shows a maximum gradient in the positive direction occurring at 
to = 13.2 rad/sec. This gradient has a value of 9a/ 33 = 0.139, while 
the corresponding o_ is _a = 0.53. In this case, the minimum singular 
value is relatively small; however, the gradient is also very small 
and therefore not very important. 

To check the validity of our conclusions, the value of 3 is 
assumed to be off by -10% (i.e., 3 = 1.8), the singular-value plot 
of the "perturbed" system is given in Figure 4.10. No significant 
change from the nominal plot (Figure 4.4) is observed. The 
stability margins have changed slightly to 
GM = -2.75 to 4.05 dB 
PM = ±21 .5° 
for o_ - 0.373. 

The results obtained above can be interpreted in the following 
manner: If the model of the plant is determined such that the 

parameter a might be expected to vary from a = _0 to a = -0.2, the 
closed-loop system may contain instabilities. Also, if the accuracy 
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FREQUENCY, to RAD/ SEC 


Figure 4.9: Singular-Value Gradient Plot that Shows the Sensitivity 

of the Relative Stability of the System with g 



Figure 4.10: Sigma Plot of the Perturbed System (a= 0 f 3 = 1.8), 

that Shows no Apparent Change of Relative Stability" 
with the Nominal System 
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on 3 is of about ±10% of the nominal value, no instabilities should 


result from such an inaccuracy. 


4,3.2 Sensitivity Analysis of the Multiloop Yaw/Roll-Pamper System 
of a Typical Business Jet 

The system analyzed in this section was discussed in Section 

3,2,3, where the controller design was given for a yaw-damper- loop 

gain of K f =1.0 and a roll-damper-loop gain of K r = 0.2» 

6 o_ 

R A 

The main objective here is to show how the sensitivity analysis 
can be used to determine the parameters that have a large effect on 
the relative .stability of MIMO closed-loop systems. 

The lateral-directional, small-perturbation equations of motion 
for the business jet are shown in Table B.l. These equations were 
used to obtain the state-space representation of the system, where 
the servo dynamics and the washout filter have been included. The 
construction of the state-space form is presented in Table 4.1, 
where the A, B, and C matrices of Equations (4.6) and (4.7) are 
shown. 

The analysis was limited to the singular-value gradients with 
respect to elements of the A matrix. Since many of the zero 
elements in the A matrix are the result of kinematic relations, they 
will not change with modelling errors (example: a(5,1) is one such 

element which is a result of the construction of the state-space 
representation of the system) ; it was therefore chosen not to, obtain 
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Table 4.1: Construction of the State-Space Representation of the Roll/Yaw Damper SAS 

for Which Sensitivity Information Was Obtained 


— — 


— 






— 


— — 


— 

— 

3 


-0.111 

0 

-0.995 

0.1594 

0.0209 

0 
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0 
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-2.133 
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0.416 
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0.3207 

0 

2.441 
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0 
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1 .0168 
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0 
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0 
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0 

0 

0 
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0 

10 
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— 



which is equivalent to 
X = AX + BU 


Equation (4.6) 




which is equivalent to 
Z - TX 

and 



which is equivalent to 


D = -CX 


Table 4.1: (Continued) 



0.2 



Table 4,1s (Continued) 


With the following representations: 



WASHOUT FILTER 



C = -0.25 £ + r 





their gradients. From the remaining elements in matrix a, only 
fourteen appeared to be important. All of these elements are 
functions of stability or control power derivatives; their gradients 
are shown in Figures 4.11 to 4.14. 

In the selection of the important parameters, several factors 
must be considered: the value of the sigma gradients, the magnitude 

of the particular element, the magnitude of the minimum singular 
value, and the value of the expected inaccuracies or changes in each 
of the elements being considered. 

To account for the magnitude of the elements and thus 
facilitate the selection, the gradient of each element was 
multiplied by the magnitude of the element. These are the * 
normalized gradients shown in Figures 4.11 to 4.14, and they can be 
interpreted as the expected change in the singular value with a 
percentage change in the element. 

Table 4.2 summarizes the information used in the final 


selection of the important elements. Column 6 gives the approximate 
percentage change in the element, p, required to obtain a 10% 
increase in the minimum singular value, a; and column 5 contains the 
normalized gradients. The difference between columns 5 and 6 is 
that column 6 takes the magnitude of the minimum singular value into 
account. This is achieved by using the following approximation: 


A 3a « 

Ap | = ~ | x — 3 — 

p p0%jj 9p/p jniax 0.10 


(4.45) 
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FREQUENCY, u (rad/sec) 


Figure 4.11: Singular-Value Gradients with Respect to Elements 

a(1,1), a(1,3), a( 1 ,4) , and a(1,5) of Matrix - of the 
Yaw/Roll Damper SAS 



Figure 4.12: Singular-Value Gradients with Respect to Elements 

a(2,1) a(2,2), a(3,1), and a (3, 5) of Matrix A of the 
Yaw/Roll Damper SAS 
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FREQUENCY, u (rad/sec) 


Figure 4.13: Singular-Value Gradients with Respect to Elements 

a(2,3), a( 2 , 5 ) , and a(2,7) of Matrix A of the Yaw/Roll 
Damper SAS 



FREQUENCY, a) (rad/secl 


Figure 4.14: Singular-Value Gradients with Respect to Elements 

a( 3 , 2) , a( 3 , 3) , and a(3,7) of Matrix A of the Yaw/Roll 
Damper SAS 



Table 4.2: Summary of the Sensitivity Analysis Information Used in the 

Selection of the Important Elements for the Yaw/Roll Damper SAS 



(1) 

(2) 

(3) 

(4) 

(5) 

(6) 


Element of A 

Frequency 

(rad/sec) 

Approximated 

a; 

3 a/ 9p 1 
— | max 

3a/ ( 3p/p) |max 

Approximated 

Ap/p |l0%o 

1) 

a( 1 , 1 ) 

0.86 

0.52 

1 

o 

. 

00 

N> 

-0.091 

-57% 

2) 

a(1 ,3) 

1.12 

0.70 

0.41 

0.41 

17% 

3) 

a(1 ,4) 

0.95 

0.57 

-1.10 

-0.17 

-34% 

4) 

a(1 ,5) 

0.72 

0.50 

-0.064 

-0.014 

-370% 

5) 

a( 2, 1 ) 

0.95 

0.57 

0.31 

0.67 

9% 

6) 

a(2,2) 

0.89 

0.53 

-0.39 

-0.21 

-25% 

7) 

a(2,3) 

0.73 

0.503 

0.05 

0.021 

240% 

8) 

a(2,5) 

0.65 

0.523 

0.022 

0.070 

75% 

9) 

a(2, 7) 

0.86 

0.52 

0.043 

0.21 

25% 

10) 

a( 3, 1 ) 

0.67 

0.52 

0.53 

0.54 

10% 

11) 

a(3,2) 

0.56 

0.55 

l 

o 

. 

in 

-0.030 

-190% 

12) 

a(3,3) ' 

1 .00 

0.60 

-0.51 

-0.095 

-65% 

13) 

a(3,5) 

0.56 

0.55 

0.059 

0.42 

1 3% 

14) 

a(3,7) 

0.49 

0.59 

0.074 

0.050 

1 20% 



The largest gradients are given by elements a(1,3) a(2,1), a(2,2), 
a.(2,7), a(3,1), a (3, 5). Although element a (1,3) has a large 
gradient, it is discarded because it is given by the following 
equation: 

Y 

a ( 1 , 3) = - 1 (4.46) 

Since >> Y r , a large change in Y r will have almost no effect on 
the value of the element; therefore, a very accurate determination 
of Y r is not warranted. This points out the importance of having a 
good understanding of the model being considered, when interpreting 
the information provided by the gradients. 

As a final check, five of the important elements (a(2,1), 
a (2, 7), a(2,2), a(3,1), and a(3,5) were allowed to change by 15% in 
their worst direction (i.e., so as to decrease jj) , and the sigma 
plot of the perturbed closed-loop system was calculated. This plot 
is shown in Figure 4.15. The minimum singular value is now given by 
a = 0. 256 with the following stability margins : 

GM = -1 .98 to 2.57 dB 
PM = ±14.7° 

Comparing the perturbed-system sigma plot (a =0.256) with the 
nominal-system sigma plot of Figure 3.25 (a = 0.502), a reduction of 
49% in the minimum singular value is observed. This reduction 
translates into reduced stability margins. The nominal system GM 
was -3.53 to 6.05 dB cbmpared to -1.98 to 2.57 dB for the perturbed 
system. Similarly, the nominal system had a PM of ±20.1° while the 
perturbed-system PM is ±14.7°. 
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Figure 4.15: Sigma Plot for System with Elements a(2,1), a(2,2), 
a{2,7), a(3,1), and a(3,5) Perturbed by 15% in Their 
Worst Direction (a = 0.256) 



Figure 4.16:' Sigma Plot for System with all A Matrix Elements 

Initially Selected Perturbed by 15% in Their Worst 
Direction (a = 0.2292) 



As a way of comparison, all fourteen elements that were 
initially selected were allowed to change by 15% (in their worst 
directions), and the corresponding sigma plot was calculated. This 
plot is shown in Figure 4.16. The minimum singular value is 
o_ ~ 0.2292 with the following stability margins: 

GM = -1 .79 to 2.26 dB 
PM = ±13.16° 

The minimum singular value has been reduced by 54% from the 
nominal system singular value, with 49% of the change due to the 
five variables selected as the most important and only 5% due to the 
remaining variables. 

A more qualitative way of looking at the same effect is given 
in the time histories of Figure 4.17 (a) and (b), where the 
difference in performance between the nominal and perturbed systems 
is noticeable; but the difference in perturbing fourteen elements, 
or only the five most important ones, is negligible. 

These results confirm the usefulness of the sensitivity 
analysis technique in the selection of the parameters that have the 
largest effect in the relative stability of the system. As has been 
mentioned before, the information obtained from the sensitivity 
analysis can be used to assign priorities for accurate determination 
of the parameters of the model. For instance, Table 4.3 shows the 
five elements selected in the previous example, along with their 
respective formulas. The stability and control power derivatives 
that affect those elements are also shown in the table. Using this 
information, a wind-tunnel test (or flight test) could be designed 
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Sideslip Angle, (degrees) Bank Angle, (degrees) 


NOMINAL SYSTEM 


FIVE ELEMENTS PERTURBED BY 152 

IN WORST DIRECTION 

— ALL ELEMENTS PERTURBED BY 152 
IN WORST DIRECTION 




Figure 4,17: Time Histories of (a) Bank Angle and (b) Sideslip 

Angle, Showing the Effect of Perturbing the Elements 
in the A Matrix of the Yaw/Roll Damper SAS by 15% in 
Their Worst Direction 


100 





Table 


Element 


a(2,1 ) 


a{ 2 , 2) 


a ( 3 , 1 ) 


a ( 3 , 5 ) 


a( 2 , 7 ) 


4.3: Summary of the Five Elements of the a Matrix 

with Largest Effect on the Relative Stability 
of the Yaw/Roll Damper SAS 


Corresponding Formula Stability or Control Power 

Derivatives 
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such that the accurate determination of C , C , C , C„ , 

n x n SL 

6 R S R 6 P 

C , C , C , and C has priority over the determination of 
l a n n. Z. 

e P 6 a 6 a 

other derivatives such as C (element 3(1,3)) and C 


(element a( 1 , 5) ) . 

In the following chapter, the results and conclusions reached 
so far are summarized. Also, some recommendations are given in the 
areas where additional research would improve the utility of the 
sensitivity analysis technique developed. 
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5. CONCLUSIONS AND RECOMMENDATIONS 


In this section, a summary of the results obtained and the 
conclusions reached is presented. Also, some suggestions are given 
regarding the kind of research activities that should be continued 
to improve the utility of the sensitivity analysis technique 
presented so far. 


5.1 SUMMARY AND CONCLUSIONS 

The primary goal of this project was the development of a 
sensitivity analysis technique for multiloop flight control 
systems. The concept of singular values, as tools in the analysis 
of relative stability of multivariable-continuous flight control 
systems, was introduced. The development of a sensitivity analysis 
technique, based on singular-value concepts, has been completed for 
continuous-multivariable flight control systems; and its use has 
been demonstrated through several examples. The application of the 
singular-value concepts to sampled-data systems has also been 
explored; a possible extension of the sensitivity analysis technique 
to cover sampled-data systems has also been proposed in Appendix C. 

A summary of general conclusions drawn from this research is 
presented below. 

1 ) The combination of the singular-value plot and eigenvalue 
plot analyses is very useful in judging the relative stability of 
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multivariable control systems. Although the stability margins 
predicted are conservative, the singular value analysis provides 
information on the potential instabilities resulting from 
simultaneous phase and gain variations in all the loops, that cannot 
always be detected by a one-loop-at-a-time frequency response 
analysis of multiloop systems. The singular value analysis provides 
the designer with the concept of multiloop stability margins, 
similar to the stability margin concept for a single-input/single- 
output system. 

2) The sensitivity analysis (^a-gradients with respect to 
system and controller parameters) technique developed is very useful 
in the identification of the parameters that have the largest impact 
on the relative stability of multiloop control systems. The direct 
relationship to changes in the stability margins of the multiloop 
system is very convenient because the control system designer has 
already developed a "feel" for the concept of phase and gain margins 
for SISO systems. 


5.2 SOME SUGGESTIONS FOR FURTHER RESEARCH 

There are five main activities which could extend and thus 
improve the utility of the sensitivity analysis technique presented 
in this report: 

1) Application of the singular-value analysis to sampled-data 
systems. The singular-value concepts and theoretical relations 
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presented are valid only for continuous systems, and they should be 
extended to handle sampled-data systems. 

2) Development of singular-value gradients for sampled-data 
systems. A suggested Approach is presented in Appendix C. 

3 ) Application of singular-value and singular-value gradient 
analyses to optimal feedback control systems, continuous and 
digital. The tools presented so far are applicable to optimal 
continuous feedback control systems. Lehtomaki, et al. [22], show 
that the stability margins for continuous-linear-quadratic (LQ) 
regulators are given by 

Gain margin, GM = 1/2, 00 
Phase margin, PM = ±60° 

and Maybeck [37] also gives the equivalent stability margins for 
digital LQ regulators. What remains to be done is to explore some 
examples of LQ regulators using the sensitivity analysis technique 
to learn how the sensitivity information could best be used. The 
concepts discussed by Cruz, et al. [16], and Safonov, et al. [17], 
on the relationship between stability and sensitivity of multiloop 
LQ regulators may prove very useful. 

4 ) Study of possible ways to reduce the conservativeness of 
the stability margins predicted by the singular values. This could 
take the form of finding another norm to define closeness to 
singularity. Another possibility is presented by Doyle [39] in the 
analysis of feedback systems with structured uncertainties. 
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APPENDIX A 


USER'S GUIDE TO SVANAL: 

THE SINGULAR VALUE GRADIENT PROGRAM 

A. 1 INTRODUCTION 

The primary purpose of SVANAL is to calculate singular value 
plots (sigma-plots) , eigenvalue plots, and singular-value gradient 
plots using the theoretical developments of Chapters Three and 
Four. The calculation of singular values is done by computing the 
singular-value decomposition (SVD) [33, 35] using the subroutine 
SNVDEC contained in ORACLS [32]. The computation of eigenvalues is 
done using the subroutine EIGEN, also contained in ORACLS [32]; and 
the computation of the singular-value gradients is done using 
Equations 4.25, 4.27, and 4.29. As an option, the user may select 
to compute the singular-value gradients with respect to elements of 
the control- law matrices using the method developed by Newsom and 
Mukhopadhyay [24] . 

The software described here consists of three main subroutines 
and several auxiliary subroutines. The flow chart in Figure A.1 
shows the relationship between* the different subroutines. The block 
named ORACLS refers to subroutines in the ORACLS package [32]. 

The user is required to write his (her) own driver program 
according to the specific needs of the task to be handled. The 
driver program must be used to read data, print or plot results, and 
set the different option flags required to calculate the sigma- 
plots, eigenvalue plots, and singular-value gradients. 
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Figure A. Is Relationship between Subroutines SVANAL, EVDIMOD, 
SVGRAD, and the Auxiliary Subroutines 
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The three main subroutines are SVANAL, EVDIMOD, and SVGRAD 


The function of SVANAL is to prepare the data according to the 
options selected by the user, calling subroutine EVDIMOD, and 
transferring input and output data from and to the user's driver 
program. 

Subroutine EVDIMOD is a modified version of a subroutine 
provided by Information and Control Systems, Inc. [36]. The three 
most important functions of EVDIMOD are the handling and 
transformation of complex matrices into augmented real matrices 
(such that algorithms and software that utilize real-matrix 
operations can be used with complex matrices); the evaluation of the 
return difference matrix (I + G 1 ) , for continuous or single-rate 
digital feedback control systems for a range of frequencies; and 
serving as the platform for the evaluation of singular values, 
eigenvalues, and singular-value gradients, through calls to ORACLS 
(SNVDEC and EIGEN) and SVGRAD subroutines, respectively. 

The main function of SVGRAD is the calculation of the singular- 
value gradients with respect to elements of the system and 
controller matrices, using the method described in Chapter Four, or 
with respect to the controller elements only, using the method 
described by Newsom and Mukhopadhyay [24]. 

In the following sections, the subroutines are described in 
more detail, the user is instructed in the use of the software, and 
an example of the kind of information that can be obtained is 
presented. 
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A. 2 SUBROUTINE EVDIMOD 

The subroutine EVDIMOD is a KU-FRL modification of the 
subroutine EIGVDI provided to us by Information and Control Systems, 
Inc. (ICS) [36]. EVDIMOD has three main functions (see also Figure 
A. 2 ) : 

1. Handling and transformation of complex matrices into 
augmented real matrices. Since this software uses subroutines from 
ORACLS [32], and ORACLS handles only real matrices, the complex 
matrices of the form (X + jY) , where X and Y are real matrices, are 
transformed into augmented real matrices of the form 



allowing us to use real matrix subroutines. 

2% Evaluation of the return difference matrix (I + ) , for 

continuous and single-rate digital feedback control systems for a 
range of frequencies. In the continuous case the plant is given by 

X = AX + BU (A. 2) 

the feedback control is 

U = -CX (A. 3) 

and the return difference matrix is 

(i + G ( s ) ] = [I + C( Is - A) ~ 1 B] (A. 4) 

where X is the state-variable vector, a and B are the system 
matrices, C is the controller matrix, U is the control vector, and I 
is the identity matrix. 
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Figure A. 2: General Flowchart of Subroutine EVOIMOD 



In the digital system case, the plant is 

v. ■ A \ * r \ <ft - 5) 

the feedback control is 

U = -CX^ (A. 6) 

and the return difference matrix is 

(i + G (z)) = [I + C( Iz - A -1 ) r] (A. 7) 

where X is the state-variable vector, U is the control vector, A and 
T are the discrete-system matrices, C is the controller matrix, and 
I is the identity matrix. 

The method used to evaluate (i + G^tz)) is presented in 
Reference [25] and is an extension of the method used to evaluate 
(l + G (s)) in the continuous case. [i + (z) ) is evaluated at 

z = e^ wT to obtain the complex matrix (i + G^ (cosu)T + j sinwT) ] , 
where j = / -1 and T is the sampling interval. EVDIMOD will 
evaluate the return difference matrix for the continuous or digital 
system, according to the option specified by the user at the time of 
data input. 

3. EVDIMOD serves as the platform for the evaluation of 
singular values, eigenvalues (of the return difference matrix), and 
singular- value gradients, through calls to ORACLS (SNVDEC and EIGEN) 
and SVGRAD . The singular values of [i + G^Ce^^)], for the digital 
system, are determined using real matrix SVD, once the complex 
matrix is transformed into the augmented real matrix form discussed 
previously. 
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Subroutine EVDIMOD uses the following auxiliary subroutines: 

INPTM: Used to input matrix A into a larger matrix, B, given a 

starting row and column number (I,J). The size of B will not 
change, 

DIAGPAR: Used to combine two matrices such that they become 

the top left and bottom right matrices of a quad partition matrix. 
The top right and bottom left matrices are set to zero in the quad 
partition matrix. 

DOUBIiM: Used to create a real (2n x 2m) matrix from a real 


(n x 2m) by doubling the number of rows: 


Real Matrix 
(n x 2m) 


(2n x 2m) 


n 

[ A1 

B1 

A2 

B2 

A3 

B3 


A1 

A2 

A3. 

. . Am 

B1 

B2 

2n 

-B1 

-B2 . , 

> . . . » 

.-Bm 

A1 

A2 


The purpose of this transformation is to do multiplication, 
addition, subtraction, and inversion of complex matrices using 
software for real matrices. Subroutines INPTM, DIAGPAR, and DOUBLM 
were provided to us by ICS. 


A. 3 SUBROUTINE SVGRAD 

The flowchart for subroutine SVGRAD is shown in Figure A. 3. 

The purpose of SVGRAD is to calculate the singular-value gradients 
with respect to elements of the system and controller matrices using 
the theoretical development of Chapter Four. As an option, the user 
can also use the method developed by Newsom and Mukhopadhyay [24] to 
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Figure A. 3: General Flowchart of Subroutine SVGRAD 
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calculate the gradients with respect to the elements of the 


controller matrices. 

For the system in Equations (A. 2) and ( A. 3 ) , SVGRAD uses the 
following expressions to evaluate the gradients with respect to 
elements in the A, B, and C matrices: 


8a. * 

= Re[$Bv.u. C<&] 

9A T 


(A. 8) 


3a. 
i 

-T 

3B 


Re[v.u. C$] 

l l 


3C T 


Re[$Bv^u^ ] 


( A.9 ) 


(A. 10) 


which correspond to Equations (4.25), (4.27), and (4.29) and 
therefore will not be discussed here again. 

It is important to note that if a discrete system has been 
selected, the gradients obtained with Equations (A. 8) and (A.9) will 
correspond to 3a^/3A T and 3a^/3T T , which are the gradients with 
respect to the matrices of the discretized system (Equation A. 5) and 
not of the original system in Equation (A. 2). Since matrices A(T) 
and T(T) are a function of matrices A and B, it is possible to find 
the gradients with respect to these latter matrices. One way in 
which this could be achieved is discussed in Appendix C; however, it 
has not been implemented at this point. 

Another option available in SVGRAD is to obtain the gradients 
using the method developed in Reference [24]. However, in order to 
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use this option some words have to be said about the plant 
(airplane) and controller representation used. 

Newsom and Mukhopadhyay [24] use the following representation 
of the plane and the feedback control law, which includes controller 
dynamics: 


X = FX + G U 
s s u 

Z = HX (A. 1 1 ) 

s 

for the plant, where *s is the state-variable vector, U is the 

control input vector, F and G u are the system matrices, Z is the 

output-variable vector, and H is the output matrix; and 

X = A X + B Z 
c c c c 

D = C X + D Z (A. 12) 

c c c 

for the controller, where X c is the controller-variable vector, and 
A c , B c , C c , and D c are the controller matrices. 

The system in Equations (A. 11) and (A. 12) can also be 
represented by Equations (A. 2) and (A. 3), if we let 


X 


A 



I 

_F_J 

B Hj 

c 1 



B = 



(A. 13) 


(A. 14) 


(A. 15) 


C 


i 

t-D H J -C ] 
c ! c 


(A. 16) 
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The expression used by Reference [24] to calculate the 


gradients is 


9 °i r 

— = Re 1 

ap A L 


- _i_ * I _ 

H( Is - A) B(v.u. )[-I i C( Is 

X 1 I 


- A) 




where p, h, and I are given by 


I 

°cj C C 

B ^ A 
C I c 

I : 

0 1 I 

I 

l 


■ [-H 


P = 


H = 


(A. 17) 


(A. 18) 


(A. 19) 


(A. 20) 


and I is the identity matrix; and and are the left and right 
normalized eigenvectors. 

The auxiliary subroutines used by SVGRAD are the following: 

DIAGPAR: Discussed in Section A. 2, 

REAX* EL: This subroutine selects the real elements of a complex 

matrix which has been stored as an augmented-real matrix in column- 
packed form: 



The returned values will correspond to matrix X in a column-packed 
form. 


DLESELM: This subroutine selects the desired elements from an 

input matrix GRAD and stores them in a matrix GRAD 2 which is the 
returned argument. 
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MAKEUV: This subroutine checks for the form in which the 


eigenvectors, u and v, are obtained from EVDIMOD and modifies them, 
if necessary, to the form required for the matrix operations in 


SVGRAD: 


input: 


output: 



or 



where the eigenvector is u = x + jy. 


A. 4 SUBROUTINE SVANAL 


The subroutine SVANAL is the executive driver of the package of 
subroutines presented in this appendix (see Figure A, 4). As such, 
it is the only subroutine which must actually be called by the 
user. SVANAL performs the following functions: 

1 . Optional echo check of all the input data 

2. Matrix and matrix dimension initialization. This includes the 
creation of the H matrix, which is used by Newsom and 
Mukhopadhyay (see Equations A. 11 -A. 19). The user inputs the H 
matrix from the equation 


Z = H X 


(A. 11 ) 


where Z is the output vector of dimension No, and Xg is the 

state vector of dimension Ns. SVNAL creates 

l 


H = 


H | 0 

V 

I 

l 


0 { I 


((NO + M) x (Ns + M) ) 
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Figure A. 4: General Flowchart* ^ u 

chart of Subroutine SVANAL 
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where H is the observer matrix defined in Equation (A.11) and I 
is an M x M identity matrix. This matrix is only created if 

A 

gradients with respect to controller parameters (p matrix) are 

A 

desired. The user can also get a printout of H if desired. 

3. Subroutine call to EVDIMOD, which handles the calculation of 
the singular values and their gradients. 

4. Output of the desired information in matrix form. The 
following information is generated by EVDIMOD (and the 
subroutines to which it is linked) : 

a) Minimum singular value and minimum eigenvalue at each 
frequency 

b) Singular value gradients with respect to selected elements 
of the system and controller matrices, at each frequency. 

The format for the input and output matrices to specify and 
retrieve the desired singular-value gradients is described in 
Section A. 5. 

A. 5 SVANAL USAGE 
A. 5.1 Calling Sequence 

CALL SVANAL (ABAR, MA, BBAR, MB, CBAR, MC, H, MH, NSM, NC, 

NO, NS, M, FREQ1 , FREQ2, NP, DELT, MI, LES, IPO, LUN1 , LUN2 , 
NWKDIM, WK, IWK, SVMAT, MSV) 
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A. 5. 2 Input Arguments 


ABAR - Variable dimensioned two-dimensional real array 

containing the augmented system matrix, A. A has 
dimensions (NSM x NSM). For discrete systems, ABAR 
must contain the state transition matrix, A. 

MA - Maximum first dimension of the array ABAR as given in 

the DIMENSION statement of the calling program. 

BBAR - Variable dimensioned two-dimensional real array 
containing the augmented control effectiveness 
matrix, B. B has dimensions (NSM x NC) . For 
discrete systems, BBAR must contain the discretized 
control power matrix, 

MB - Maximum first dimension of the array BBAR as given in 

the DIMENSION statement of the calling program. 

CBAR - Variable dimensioned two-dimensional real array 
containing the augmented feedback law matrix, C. 

C has dimensions (NC x NSM). 

MC - Maximum first dimension of the array CBAR as given in 

the DIMENSION statement of the calling program. 

H - Variable dimensioned two-dimensional real array 

containing the augmented observer matrix, H. This 
matrix is not required unless gradients with respect 
to parametrs in P (Newsom and Mukhopadhyay method) 
are desired, but the variable name must appear in the 
call statement. 
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MH 


NSM 

NC 

NO 

NS 

M 


NOTE: 


FREQ1 

FREQ2 

NP 


DELT 


MI 


Maximum first dimension of the array H as given in 
the DIMENSION statement of the calling program 
Dimension of the augmented state-variable vector. 
Dimension of the control vector. 

Dimension of the output vector. 

State-vector dimension before augmentation. 

Order of the controller-dynamics equation. See 
Equation (A. 12). 

The variables MH, NO, NS, and M must appear in the 
call statement to SVANAL; but they do not have to be 
initialized if gradients with respect to parametrs in 

A 

p are not desired. NO, NS, M, and NC are consistent 
with the definitions given in Reference [23] and 
Section A. 3. 

Lowest frequency for the singular value analysis. 
Highest frequency for the singular value analysis. 
Number of frequencies between FREQ1 and FREQ2 where 
singular values, eigenvalues, and gradients are to be 
calculated. 

Sampling interval for discrete system. This variable 
must be set to zero for continuous systems. 

Integer vector whose dimension is 4. Each entry 
indicates the number of gradients desired from the 
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A, B, C, and P matrix, respectively. This input 
argument is described in more detail in Section 

A. 5 .8. 

LES - Integer array whose first dimension must be 4. 

Second dimension must be max k (MI(k)). Entry (k,i) 
indicates the location (counting across the rows) in 

A 

the k-th matrix (a, B, C, or P) of the parameter with 
respect to which gradient plot is desired. This 
input argument is described in more detail in Section 

A. 5. 8. 

IPO - Integer variable indicating the type of output 
desired. 

IPO = 0: All printout suppressed. 

A 

IPO = 1 : Input matrix data echo; includes H matrix 

A 

printout if P gradients are requested. 

This printout is directed to logical unit 
number LUN1 • 

IPO = 2: Data file printout. This printout is 

directed to logical unit number LUN2 and is 
described in Section A. 5. 8. 

IPO = 3: Both matrix data echo and data file 

information are printed out to their 
respective logical unit numbers. 

LUN1 - Logical unit number for the matrix data echo and 

error output. This variable should be set to the LUN 
for the terminal. 
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LUN2 - Logical unit number for the data file printout (see 
also Section A, 5. 8). The calling program must open 
and initialize the data file referenced by LUN2. 

NOTE: LUN1 must be initialized to allow error statements. 

LUN2 need not be initialized if IPO < 2. 

NWKDIM - Dimension of the work vector, WK, as given by the 

DIMENSION statement of the calling program. NWKDIM 
must be at least 

2 2 
NWKDIM = 37NSM + 36NSM + 6(NSM • NC) + 28NC + 42NC + 93 + 

+ 5MAX ( 4NSM , 4(N0 + M) (NC + M) ) + MAx(MI(k), k = 1,4) + 

+ 4 (NO + M) (NSM + NC) + 4(NSM • M) (A. 21) 

WK - Real work vector dimensioned at least NWKDIM in the 
calling routine. 

IWK - Integer work vector dimensioned at least 2NC + 2M in 
the calling routine. 

MSV - Maximum first dimension of the array SVMAT as given 
in the DIMENSION statement of the calling program. 

MSV must be greater than NP. 
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A. 5 . 3 Output Arguments 


SVMAT - Two-dimensional array containing all output 

information. SVMAT must be dimensioned at least 
[NP x (Mid) + MI (2) + MI ( 3 ) + MI (4) + 3j] 
in the calling routine. Data is organized as 
follows: 

SVMAT (I, 1): frequency, in rad/sec, 1=1 to NP 

SVMAT (I, 2): minimum singular value, 1=1 to NP 

SVMAT (I, 3): minimum eigenvalue, 1=1 to NP 

SVMAT (I, J + 3): gradients with respect to 

parameters in the A matrix, J = 1 to- Mid ) 

1 = 1 to NP 

SVMAT (I, J -t MI(1) + MI ( 2) + 3): gradients with 

respect to parameters in the C matrix, J = 1 to MI(3) 

1 = 1 to NP 

SVMAT (I, J + MI ( 1 ) + MI ( 2) + Ml(3) + 3): 

A 

gradients with respect to parameters in the p matrix, 
J = 1 to MI (4) 

1=1 to NP 


128 



A. 5*4 COMMON Blocks 


The COMMON blocks employed by SVANAL are 
GRADSTF and SAVE, 


A. 5.5 Error Messages 


(1) If the work array is not large enough, the following message 

will be printed: "THE WORK ARRAY IS NOT LARGE ENOUGH IN 

SNGVDI." 

(2) If errors occur within SNVDEC, the following message will be 

printed: "PASS # THERE IS AN ERROR IN THE SECOND CALL 

TO SNVDEC, I ERR = IN SNGVDI." 


A. 5.6 Subroutines Employed by SVANAL 


The following subroutines from ORACLS, ICS, and FRL are 
employed : 

ORACLS - ADD, NULL, MULT, JUXTR, JUXTC, UNITY, EQUATE, PRNT, 
MAXEL, AXPXB, SCALE, SNVDEC, HSHLDR, BCKMLT, SCHUR, 
SHRSLV , SYSSLV, NORMS, LNCNT, EIGEN, INVIT, BALBAK, 
ELMBAK, HQR, ELMHES, BALANC 
ICS - DOUBLM, DIAGPAR, INPTM 

FRL - EVDIMOD, SVGRAD, DLESELM, MAKEUV, REALEL 
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A. 5. 7 Subroutines Employing SVANAL 


None • 


A. 5.8 Concluding Remarks 

A singular value gradient exists with respect to each element 

A 

of the A, B, C, and P matrix (these matrices are described in 
Section A* 3). SVANAL outputs only those gradients which are 
specified by the user. To identify the gradients of interest, the 
following information must be given to SVANAL: 

1 ) MI(k) f k = 1 , 2, 3, and 4. 

MI is a vector containing four integers. Each integer 
represents a particular matrix: 


mi(i ) 

====> 

A matrix 

MI( 2) 

”> 

B matrix 

MI(3) 

— > 

C matrix 

MI(4) 

====> 

P matrix 


The value of MI(k) indicates the number of gradients 
desired from the corresponding matrix. For instance, 
MI (3) =7 indicates that 7 gradients with respect to 
variables in the c matrix are desired. 

2) LES (k , j ) , k = 1, 2, 3, and 4 

j = 1, 2, 3,..., MI (k ) 
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LES is an integer matrix containing four rows and enough 


columns to accommodate the largest value of MI(k). Each 
row corresponds to one of the four matrices, 

/■v 

A, B, C, and P, respectively (in the same order as the 

elements of MI(k)). Each row contains exactly MI(k) 

nonzero elements, one for each of the desired gradients. 

The value of each element indentifies the matrix location 

(in the corresponding A, B, C, or P matrix) of the 

parameter for which a gradient plot is desired. Matrix 

locations are numbered from left to right across the rows, 

continuing on the next row after each row is numbered. 

For instance, a 4 x 4 matrix is numbered as follows: 

12 3 4 

5 6 7 8 

9 10 11 12 

13 14 15 16 

As an example, suppose the A matrix is 4 x 4. Note that the 
dimension of the array ABAR, in which A is stored, is of no 
consequence. If the user wants gradients with respect to elements 
2, 7, 9, and 12 (as numbered above), he would set Ml(1) = 4 
(indicating that four gradients are desired from the A matrix) and 
LES ( 1 , 1 ) = 2 
LES ( 1 , 2) = 7 
LES ( 1 , 3) = 9 
LES ( 1 , 4) =12. 
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Refer to Section A. 6 for a complete example of how MI and LES are 
created. 

Data-file output is provided to allow part of the information 
generated by SVANAL to be saved in the instance of a program 
abort. Numerical problems in the ORACLS subroutine SNVDEC sometimes 
cause errors, at the time of execution, which cannot be avoided. If 
the minimum singular value, minimum eigenvalue, and gradient are 
printed out at each frequency as they are calculated (IPO = 2 or 3), 
then this information will not be lost if an error occurs. The 
plots can then be continued at or near the frequency for which the 
error occurred. The user then proceeds, segment by segment, across 
the frequency range until the entire plot has been obtained. 

The format of the output file is best described by the WRITE 

and FORMAT statements which create it. At each frequency, FREQ, 

SVMIN, and EVMIN are the minimum singular value and minimum 

eigenvalue, respectively. LGRAD( I, K) contains the gradients, where 

/\ 

1=1 to 4 represent each of the matrices A, B, C, and P in that 

order. K goes from 1 to MI (I), the number of gradients desired from 

the respective matrix. 

WRITE ( LUN2 , 4000 ) FREQ, SVMIN, EVMIN 
DO 1 I =1 ,4 

1 WRITE (LUN 2,4000) ( LGRAD ( I , K ) K = 1 , MI ( I ) ) 

4000 FORMAT (8F10.5) 
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A, 6 EXAMPLE 


In Section 4.3*. 1, a singular value analysis was done on the 
following single-input/single-output system. Two cases were 
examined, one in which the singular-value gradient with respect to a 
parameter in the C matrix was calculated, and one in which gradients 
with respect to parameters in the A matrix were calculated. In both 
of these cases, however, the nominal system is the same: 




— 


- 


- — 


- - 

XI 


0 

1 

0 


xl 


0 

. 









x2 

- 

0 

0 

1 


x2 


0 

x3 


-40 

-28 

-6 


x3 


1 





. 


— 


- - 



In case I, the gradient with respect to element 1 in C was 
obtained. In case II, gradients with respect to elements 7, 8, and 
9 in A were calculated. The frequency range for all the plots in 
Section 4.3.1 was .01 to 100. rad/sec. All four of these gradients, 
as well as the minimum singular value and minimum eigenvalue plots, 
can be found using one single call to SVANAL, with the following 
input information: 
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ABAR = 


0 
0 

L-40 


1 0 
0 1 
-28 -6. 


BBAR = 


0 

0 

Li. 


CBAR = £ 0 200 oj 

MA, MB, MC: Depend on dimension statement of calling program. 

A 

H, MH: Not needed; p gradients not being calculated. 

NSM = 3 


NC = 1 

NO, NS, M: Not needed. 

FREQ1 = .01 
FREQ2 = 100. 

NP = 150 

DELT =0.0 (Continuous System) 

MI(1) = 3 
MI( 2) = 0 
MI( 3) = 1 
MI( 4) = 0 

LES ( 1 , 1) = 7 
LES ( 1 , 2) = 8 
LES( 1 , 3) = 9 
LES( 3, 1) = 1 

IPO = 3 

LUN1 = 3 (LUN for the CRT) 

LUN2 = 18 (LUN for the data file) 

NWKDIM = 700 minimum required according to Equation (A. 21) 
MSV: Depends on dimension statement of calling program. 
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The output information all appears in the matrix SVMAT. 
Table A.1 is an image of this matrix, with the contents of each 
colunn labeled to indicate how output information is arranged in 
SVMAT, 
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Table A.1: An Example of Output-Matrix SVMAT 

Containing Information on Frequency, 
Minimum Singular- and Eigenvalues, 
and Singular-Value Gradients 


SVMAT ( I , J ) : 


NP 


1 

o 

3 

4 

5 

6 

7 

□mega 

MinSig 

MinEig 

dSig/dFl 

dSig/dF2 

dSi g/dF3 

dSig/dAlf 

0. 10000 

1.14864 

1.14864 

0.00693 

-0.00104 

-0.00007 

0.02169 

0. 10965 

1. 17636 

1 . 1 7636 

0.00812 

-0.00121 

-0.00010 

0.02117 

O. 12023 

1.20881 

1.20881 

0.00949 

-0.00139 

-0.00014 

0.02058 

0. 13183 

1.24667 

1 . 24667 

0.01105 

-0.00160 

-0.00019 

0.01994 

0. 14454 

1.29067 

1 . 29067 

0.01281 

-0.00183 

-0.00027 

0.01923 

0. 1SB49 

1.34158 

1.34158 

0.01478 

-0.00207 

-0.00037 

0.01848 

0. 17378 

1.40023 

1 . 40023 

0.01698 

-0.00232 

-0.00051 

0.01767 

0. 19055 

1.46751 

1.46751 

0.01941 

-0.00258 

-0.00070 

0.01683 

0.20893 

1.54432 

1.54432 

0.02209 

-0.00233 

-0.00096 

0.01595 

0.22909 

1.63162 

1.63162 

0.02502 

-0.00306 

-0.00131 

0.01505 

0.251 19 

1 . 73039 

1 . 73039 

0.02820 

-0.00325 

-0.00178 

0.01414 

0.27542 

1.84166 

1.84166 

0.03164 

-0.00338 

-0.00240 

0.01323 

0.30200 

1 . 96646 

1 . 96646 

0.03534 

-0.00342 

-0.00322 

0.01233 

0.33113 

2. 10587 

2. 10587 

0.03929 

-0.00332 

-0.00431 

0.01144 

0.36308 

2.26097 

2.26097 

0.04349 

-0.00304 

-0.00573 

0.01058 

0.3981 1 

2.43286 

2. 432B6 

0.04792 

-0.00249 

-0.00759 

0.00974 

0.43652 

2. 62264 

2.62264 

0.05255 

-0.00160 

-0.01001 

0.00894 

0.47863 

2.83139 

2.B3139 

0.05737 

-0.00025 

-0.01314 

0.00818 

0.52481 

3.06017 

3.06017 

0.06231 

0.00171 

-0.01716 

0.00745 

0.57544 

3.30996 

3.30996 

0.06732 

0.00444 

-0.02229 

0.00676 

0.63096 

3.58163 

3.58163 

0.07230 

0.00818 

-0.02878 

0.00611 

0.69183 

3.B7593 

3.87593 

0.07715 

0.01319 

-0.03693 

0.00550 

0.75SSB 

4. 19341 

4. 19341 

0.08173 

0.01977 

-0.04703 

0.00493 

0.83176 

4.53439 

4.53439 

0.08586 

0.02831 

-0. 05940 

0.00439 

0.91201 

4.89892 

4.89B92 

0.08935 

0.03921 

-0.07432 

0. 003B8 

1.00000 

5.28674 

5.28674 

0.09195 

0.05295 

-0.09195 

0.00341 

- 

- 

- 

• 

- 

- 

• 

• 

• 

• 

• 

• 

• 

• 
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APPENDIX B 


AERODYNAMIC AND GEOMETRIC DATA FOR A 
TYPICAL BUSINESS JET 

This appendix contains the geometric and lateral-directional 
aerodynamic data used in the example analyses of the yaw-damper and 
yaw-roll-damper SAS of a typical business jet. The data were 
obtained from Reference [31], 

The lateral-directional, small-perturbation equations of motion 
of the airplane are shown in Table B.l. The formulas for the 
lateral-directional dimensional stability derivatives were also 
obtained from Reference [31] and are contained in Table B.2. The 
dimensionless derivatives and the initial characteristics of the 
airplane for approach flight condition and full-flaps down at sea. 
level are shown in Table B.3* 
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Table B.1: Lateral-Directional, Small-Perturbation Equations 

of Motion for a Business Jet Airplane 




Table B.2: Lateral-Directional Dimensional Stability Derivatives 

(from Reference [31]) 
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Table B,3: Aerodynamic and Inertial Characteristics Needed 

for the Yaw- and Yaw/Roll -Damper Analyses 
for a Typical Business Jet (from Reference [31]) 


The following data are for a business jet airplane in an approach 
flight condition. 


Full Flaps - Sea level 


Vg = 120 kts 

U l = 

202 fps 


= 23,500 slug-ft 2 

W ® 11,800 lbs 




s 

S = 232. ft 2 




[ = -100 slug-ft 2 

xz ° 





s 

b = 34.2 ft 




[ = 40,700 slug-ft 2 

zz 

s 

C L = 1.05 

M = 

.18 


C = - 0.73 rad" 1 

C 

= + 0.4 rad 1 

C =0 

-1 

y r 



y p 

Co = - 0.13 rad 

s 

C * 

=0.30 rad 

= -.386 rad 

r 



P 

C = 0.107 rad" 1 

n B 

C 

n 

r 

= -0.233 

rad 1 

C = - 0.066 rad" 1 
n 

P 

C = .1375 rad" 1 


C 

= 0 


y 6 


y s 



R 


A 



C = .0194 rad” 1 

X'x 


% 

= 0.149 rad 1 

6 r 


6 a 



'c = -.0745 rad" 1 


C 

= -.0355 rad’ 1 

n 6 


n 6 



R 


A 




Yaw-Rate-To- Rudder Transfer Function: 


I_ (s) = - 0 708 (s + 0.912) [(8 - .136) 2 + (.521) 2 ] 

X ' i , , n OTOW- _ A A A \ r f ~ l A A O 1 \ 9 l / 1 


(s + 0.832) (s - .004) [(s + .0031)2 + (1.10)2] 


Roll-Rate-To-Aileron Transfer Function: 


s[(s + .120) 2 + (0.846) 2 ! 


o ^ a 2 44 

6 a K ‘ (s + 0.832) (s - .004) [(s + .0031)2 + ( 1 . 10 ) 2 ] 
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APPENDIX C 


A SUGGESTED APPROACH TO THE DEVELOPMENT 
OF SINGULAR-VALUE GRADIENTS 
FOR SAMPLED-DATA SYSTEMS 

This appendix presents a suggested approach to the development 
of the singular-value gradients for sampled-data systems. The use 
of these expressions should extend the applicability of the 
sensitivity analysis technique to sampled-data systems. 

Sampled-data systems are composed of a continuous dynamical 
subsystem and a discrete subsystem. The continuous subsystem is 
called the plant, and the discrete subsystem is called the digital 
controller. To analyze a sampled-data system, the continuous plant 
must be discretized so that the two subsystems have a common 
representation. If the continuous model of the plant is given by 
Equation (4.6), the following equation represents the discretized 
model : 

x k+i = a(t)x^ + r(T)u k (C.1) 

where X is the state-variable vector, U is the control vector, and 

A(T) is called the transition matrix, given by 

2 

A(T) = I + AT + A 2 y- + . . . (C.2) 

and T(T) is given by 
T 

r(T) = [ /A(t)<3t]B (C. 3) 

O 

where A and B are the continuous system matrices, I is the identity 
matrix, and T is the sampling period. The digital feedback control 
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law is given by 


U = -C^ (C. 4) 

Therefore, if the expressions for singular-value gradients of 
Section 4.1 are applied, the gradients obtained will be with respect 
to elements of the discretized matrices, A(T) and T(T). In order to 
evaluate the gradients with respect to the original matrices A and 
B, a different set of expressions must be developed. These 
expressions can be developed following the same technique used in 
Section 4.1; but since A(T) and T(T) are expressed in terms of 
series, some complications arise and a decision must be made as to 
how many terms of the series must be carried. The following is a 
description of the suggested approach, where only first- and second- 
order terms have been considered in Equations (C.2) and (C.3). In 
this case, Equation (C.3) can be expressed as 

_ t 2 -2 T^ - 

r(T) = [IT + A + A “ — ] B ( C. 5 ) 

2 6 

The general equation of the singular-value gradients (Equation 
4.17) can now be written as 


3o. 

i 

3p 


Re 


, r 3d + c*r> 

tr [ 

3p 


★ 

v.u. ] 

1 X 


(C.6) 


where $ is now given by 


$ = (Iz - A) 


ioj. T 
i 

z=e 


Expanding Equation (C.6), we obtain 


(C. 7 ) 
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• J}p — £ r» __ A * 

— — = Re • tr[(— OF + CO — + CO — ODv.u. ] 
dp dp dp dp 1 1 


(C.8 ) 


and it is possible to obtain three expressions for do, /dp, one each 
for the following cases: 

/S A 

a) Let p be an element of p and p = B; then we obtain 


9a i - 3T * 

t — = Re • tr[(C$ — )v.u. ] 
9p 3p i 1 


(C. 9) 


with 


— -2 ~ A 
T = [IT + A — — + A — ]P = A'P 
2 6 


(C.10) 


Equation (C.9) is then written as 


3a. _ a ; * 

-r-i = Re • tr [ (C$A' ) -r— (v.U. )] 

3p dp X X 


(C.11 ) 


and, by using Equations (4.22) and (4.23), we obtain 


_ * * * 
= Re ( (C$A' ) (v.u. ) ] 

3P 11 


(C. 1 2) 


or, using the transpose. 


3a. 3a. * 

~ - — ~ = Re[ (v.u. )(C*A')1 (C. 1 3) 

3P 3b 1 1 


b) Let p be an element of p and p = a, then we obtain from 


Equation (C.8) 


da. 

i 

3 P 


Re 


— 3 P 

tr[(C$ ir + 

3p 


C4 -jp- 4-Dv.u.*] 
dp 11 


with 


(C. 1 4 ) 


2 3 

rri A rn J ^ A _ 

r = [it + — p + — pp]b 

2 6 
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A=I+TP+— PP 

Substituting (C.15) in (C.14), and rearranging, we obtain 


9a. 


„2 „ 


„3 A . 


8P 9P 


= — Re • tr[{(C$ — )P(B) + (C$ g-)PP(B) + 


+ (C*T)P(*D + (C$ |-)pp($r)}(v i u i )] 


Using Equation (4.23) and 


(C.15) 


(C. 16) 


3 T T T T T T 

— trace [BAAC] = B C A + A B C 
3A 


(C.17) 


we obtain, using the transpose 


3a. 3a. _ 2 __ 3 

~~ = — - = Re [ ( Bv.u. CO — ) + (ABv.u. CO — ) + 
* T -T x x 2 x x 6 

3P 3a 


*- t _ *_ 

+ (Bv.u. CO — A) + (OTv.u. COT) + 
x x 6 XX 


c) 


2 ^ 2 

- + (AOTv.u. CO ~~) + (OTv.u. CO — A)] (C.18) 

x x 2 x x 2 

Let p be an element of p and p = C; then, from Equation 


( C. 8 ) 

^ a i 3C * 

^ = Re • tr[—- OTv.u. ] ( C * 1 9 ) 

3p 3p x x 

and 


3a. 3a. 

i x 



* 

RefOTv^u^ ] 


(C. 20) 


which is similar to Equation ( 4 • 29 ) with T instead of B, and the 
discrete version of 0. 
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